函数文件1：real_fun.m
1 function f=real_fun(x0,t0) 2 %精确解 3 f=4*x0*(1-x0)*sin(t0);
函数文件2：F.m
1 function f=F(N,u,U,t,h1,h2) 2 %非线性方程组 3 %h1是x的步长，h2是t的步长 4 %u表示迭代节点,上一时刻的数值解 5 %h表示时间节点上的步长 6 %N表示空间节点的步数 7 a0=0.5*t^4*h2*N^2; 8 f(1,1)=a0*(U(2)^2-2*U(1)^2)+h2*fi(h1,t)+u(1)-U(1); 9 f(N-1,1)=a0*(-2*U(N-1)^2+U(N-2)^2)+h2*fi((N-1)*h1,t)+u(N-1)-U(N-1); 10 for p=2:N-2 11 f(p,1)=a0*(U(p+1)^2-2*U(p)^2+U(p-1)^2)+h2*fi(p*h1,t)+u(p)-U(p); 12 end
函数文件3：fi.m
1 function f=fi(x0,t0) 2 %等式右边的f函数 3 f=4*x0*(1-x0)*cos(t0)-16*t0^4*(6*x0^2-6*x0+1)*(sin(t0))^2;
函数文件4：Jacobian.m
1 function g=Jacobian(n,u,t,h1,h2) 2 %计算每个时间节点的牛顿迭代过程中的雅可比矩阵 3 %u表示迭代初值，上一时刻的数值解作为迭代初值 4 a=0.5*t^4*h2*n^2; 5 g=zeros(n-1); 6 g(1,2)=2*a*u(2); 7 g(1,1)=-4*a*u(1); 8 g(n-1,n-1)=-4*a*u(n-1); 9 g(n-1,n-2)=2*a*u(n-2); 10 for p=2:n-2 11 g(p,p+1)=2*a*u(p+1); 12 g(p,p)=-4*a*u(p); 13 g(p,p-1)=2*a*u(p-1); 14 end 15 g=g-eye(n-1);
函数文件5：Newtond.m
1 function x=Newtond(n,u,t,h1,h2) 2 %使用修改后的牛顿迭代，可以不求雅可比de逆 3 %U中间代初值 4 %u起始迭代初值 5 U=u; 6 tol=0.5e-5; 7 % Jacobi=Jacobian(n,u,t,h1,h2);%每隔k步求一次雅可比 8 x1=U-Jacobian(n,u,t,h1,h2)\F(n,u,U,t,h1,h2); 9 while (norm(x1-U,1)>=tol) 10 %数值解的1范数是否在误差范围内 11 U=x1; 12 x1=U-Jacobian(n,u,t,h1,h2)\F(n,u,U,t,h1,h2); 13 end 14 x=x1;%不动点
脚本文件：
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 ti=0:t_h:0.5;%t的节点 10 %********************真解************************** 11 for i=1:length(x) 12 for j=1:length(ti) 13 real_Z(i,j)=real_fun(x(i),ti(j)); 14 end 15 end 16 %********************真解************************** 17 %********************数值解************************** 18 ui=zeros(length(x)-2,1);%牛顿迭代初值 19 Z=zeros(length(x),length(ti)); 20 for i=1:length(ti)-1 21 Z(2:length(x)-1,i+1)=Newtond(length(x)-1,ui,ti(i+1),x_h,t_h);%t(i+1)时间的牛顿数值解 22 ui=Z(2:length(x)-1,i+1);%牛顿迭代初值，上一时刻的数值解作为迭代初值 23 end 24 25 %********************数值解************************** 26 [X,Y]=meshgrid(x,ti); 27 subplot(2,2,1), 28 mesh(X,Y,real_Z'); 29 xlabel('x');ylabel('t');zlabel('u');title('analytical solution'); 30 subplot(2,2,2), 31 mesh(X,Y,Z'); 32 xlabel('x');ylabel('t');zlabel('u');title('numerical solution'); 33 subplot(2,2,3), 34 mesh(X,Y,real_Z'-Z'); 35 xlabel('x');ylabel('t');zlabel('u');title('error solution'); 36 title('牛顿迭代法'); 37 grid on; 38 toc;
效果图：
