Неоднородное уравнение теплопроводности - PullRequest
Новая
       17

Неоднородное уравнение теплопроводности

0 голосов
/ 30 апреля 2019

У меня есть этот код, который решает уравнение теплопроводности с неоднородными граничными условиями: $$ u_t = ku_ {xx} $$ где k - коэффициент диффузии но теперь я хочу добавить константу, которая делает уравнение неоднородным, например так: $$ u_t = ku_ {хх} + H $$ 

function [u,er]=ftcs(t0,tf,nt,a,b,nx,ci,cca,ccb,dif,cada,sol)
% explicit method

%u solution
%er vector error
%t0 y tf time limits
%nt 
%nx 
%a b x value
%ci initial condition
%cca y ccb boundary conditions
%dif diffusion coefficient
%cada time lapse
%sol solution if it's known

x=linspace(a,b,nx);  x=x';  dx=x(2)-x(1); %h
t=linspace(t0,tf,nt); t=t';  dt=t(2)-t(1); %k
r=dt/dx^2*dif;  
if r>.5
    clc
    disp('stability not satisfied')
    pause
end
di=1-2*r;
cca=inline(cca,'t');
ccb=inline(ccb,'t');
A=diag(di*ones(nx-2,1))+diag(r*ones(nx-3,1),1)+diag(r*ones(nx-3,1),-1);
A=[[r;zeros(nx-3,1)] A [zeros(nx-3,1);r]];
ci=vectorize(inline(ci,'x')); 
u=ci(x);      % vectorwith a solution for t=0, initial condition
gra=plot(x,u);
axis([a b min(u) max(u)])
pause
z=u;
u(1)=cca(t(1)); 
u(end)=ccb(t(1)); 
for i=1:nt-1;  
    u=A*u;
    u=[cca(t(i+1)); u ; ccb(t(i+1))];
    if mod(i,cada)==0
        set(gra,'ydata',u');
        pause(.1)
        z(:,end+1)=u;
    end
end
if nargin==12 & nargout==2
    sol=vectorize(inline(sol,'x','t'));
    er=u-sol(x,tf);
end

pause
figure
surfc(z);shading interp; colormap(hot);set(gca,'ydir','reverse');
rotate3d
...