Matlab:线性热传导(抛物线方程)问题

Matlab:线性热传导(抛物线方程)问题

Matlab:线性热传导(抛物线方程)问题

函数文件1:real_fun.m

 function f=real_fun(x0,t0)
f=(x0-x0^2)*exp(-t0);

函数文件2:fun.m

 function f=fun(x0,t0)
f=(x0^2-x0)*exp(-t0)+2*exp(-t0);

函数文件3:fi.m

 function f=fi(x0)
f=x0-x0^2;

脚本文件:

 tic;
clc
clear
N=100;
M=1000;
t_h=1/M;%t的步长
x_h=1/N;%x的步长
x=0:x_h:1;%x的节点
t=0:t_h:1;%t的节点
B=-2*ones(1,N-1);
C=1*ones(1,N-2);
D=1*ones(1,N-2);
A=diag(B)+diag(C,1)+diag(D,-1);%三对角矩阵
F=zeros(N-1,M);
for i=1:N-1
for j=1:M
F(i,j)=fun(x(i+1),t(j));
end
end
F=F.*t_h;
%********************数值解************************************
J=-N^2*A*t_h+eye(N-1);%求解线性方程组的系数矩阵
initial=zeros(N-1,1);
z=zeros(N-1,M); for i=1:N-1
initial(i)=fi(x(i+1));
end
b=initial;
for j=1:M%控制t的节点
a=b;
a=a+F(1:N-1,j);
z(1:N-1,j)=J\a;%解是n-1维的
b=z(1:N-1,j);%变成下一次求解的初值
end
z=[initial,z];
Z=[zeros(1,M+1);z;zeros(1,M+1)]; %********************数值解************************************
for i=1:N+1
for j=1:M+1
real_Z(i,j)=real_fun(x(i),t(j));
end
end
compare=abs(real_Z-Z);
[X,Y]=meshgrid(x,t);
% colormap(jet) subplot(2,2,1),
mesh(X,Y,Z');
xlabel('x');ylabel('t');zlabel('u');title('analytical solution');
subplot(2,2,2),
mesh(X,Y,real_Z');
xlabel('x');ylabel('t');zlabel('u');title('numerical solution');
subplot(2,2,3),
mesh(X,Y,compare');
xlabel('x');ylabel('t');zlabel('u');title('error solution');
grid on;
toc;

效果图:

Matlab:线性热传导(抛物线方程)问题

上一篇:Java获取主板序列号、MAC地址、CPU序列号工具类


下一篇:OC基础--构造方法 id类型