1 function f=real_fun(x0,t0)
2 %精确解
3 f=4*x0*(1-x0)*sin(t0);

 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

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;

 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);

 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;