函数文件1：real_fun.m
1 function f=real_fun(x0,t0) 2 f=(x0-x0^2)*exp(-t0);
函数文件2：fun.m
1 function f=fun(x0,t0) 2 f=(x0^2-x0)*exp(-t0)+2*exp(-t0);
函数文件3：fi.m
1 function f=fi(x0) 2 f=x0-x0^2;
脚本文件：
1 tic; 2 clc 3 clear 4 N=100; 5 M=1000; 6 t_h=1/M;%t的步长 7 x_h=1/N;%x的步长 8 x=0:x_h:1;%x的节点 9 t=0:t_h:1;%t的节点 10 B=-2*ones(1,N-1); 11 C=1*ones(1,N-2); 12 D=1*ones(1,N-2); 13 A=diag(B)+diag(C,1)+diag(D,-1);%三对角矩阵 14 F=zeros(N-1,M); 15 for i=1:N-1 16 for j=1:M 17 F(i,j)=fun(x(i+1),t(j)); 18 end 19 end 20 F=F.*t_h; 21 %********************数值解************************************ 22 J=-N^2*A*t_h+eye(N-1);%求解线性方程组的系数矩阵 23 initial=zeros(N-1,1); 24 z=zeros(N-1,M); 25 26 for i=1:N-1 27 initial(i)=fi(x(i+1)); 28 end 29 b=initial; 30 for j=1:M%控制t的节点 31 a=b; 32 a=a+F(1:N-1,j); 33 z(1:N-1,j)=Ja;%解是n-1维的 34 b=z(1:N-1,j);%变成下一次求解的初值 35 end 36 z=[initial,z]; 37 Z=[zeros(1,M+1);z;zeros(1,M+1)]; 38 39 %********************数值解************************************ 40 for i=1:N+1 41 for j=1:M+1 42 real_Z(i,j)=real_fun(x(i),t(j)); 43 end 44 end 45 compare=abs(real_Z-Z); 46 [X,Y]=meshgrid(x,t); 47 % colormap(jet) 48 49 subplot(2,2,1), 50 mesh(X,Y,Z'); 51 xlabel('x');ylabel('t');zlabel('u');title('analytical solution'); 52 subplot(2,2,2), 53 mesh(X,Y,real_Z'); 54 xlabel('x');ylabel('t');zlabel('u');title('numerical solution'); 55 subplot(2,2,3), 56 mesh(X,Y,compare'); 57 xlabel('x');ylabel('t');zlabel('u');title('error solution'); 58 grid on; 59 toc;
效果图：
