Matlab:非线性高阶常微分方程的几种解法
一、隐式Euler：
函数文件1：
1 function b=F(t,x0,u,h) 2 b(1,1)=x0(1)-h*x0(2)-u(1); 3 b(2,1)=x0(2)+2*h*x0(2)/t+4*h*(2*exp(x0(1))+exp(x0(1)/2))-u(2);
函数文件2：
1 function g=Jacobian(x0,t,h) 2 g(1,1)=1; 3 g(1,2)=-h; 4 g(2,1)=4*h*(2*exp(x0(1))+0.5*exp(x0(1)/2)); 5 g(2,2)=1+2*h/t;
函数文件3：
1 function x=newton_Iterative_method(t,u,h) 2 % u：上一节点的数值解或者初值 3 % x0 每次迭代的上一节点的数值 4 % x1 每次的迭代数值 5 % tol 允许误差 6 % f 右端函数 7 x0=u; 8 tol=1e-14; 9 x1=x0-Jacobian(x0,t,h)\F(t,x0,u,h); 10 while (norm(x1-x0,inf)>tol) 11 %数值解的2范数是否在误差范围内 12 x0=x1; 13 x1=x0-Jacobian(x0,t,h)\F(t,x0,u,h); 14 end 15 x=x1;%不动点
脚本文件：
1 tic; 2 clear all 3 clc 4 N=[128 64 32 16 8 4]; 5 for j=1:length(N) 6 h=1/N(j); 7 x=0:h:1; 8 y=inline('-2*log(1+x.^2)','x'); 9 Accurate=y(x); 10 Numerical=zeros(2,N(j)+1); 11 for i=1:N(j) 12 Numerical(1:2,i+1)=newton_Iterative_method(x(i+1),Numerical(1:2,i),h); 13 end 14 error(1:N(j)+1,j)=Numerical(1,:)-Accurate; 15 figure(j) 16 subplot(1,3,1) 17 plot(x,Accurate); 18 xlabel('x'); 19 ylabel('Accurate'); 20 grid on 21 subplot(1,3,2) 22 plot(x,Numerical(1,:)); 23 xlabel('x'); 24 ylabel('Numerical'); 25 grid on 26 subplot(1,3,3) 27 plot(x,error(1:N(j)+1,j)); 28 xlabel('x'); 29 ylabel('error'); 30 title(1/N(j)); 31 grid on 32 end 33 for k=2:length(N) 34 X=norm(error(:,k),inf)/norm(error(:,1),inf); 35 H=N(1)/N(k); 36 Y(k-1)=log(X)/log(H); 37 end 38 figure(length(N)+1) 39 plot(1:length(N)-1,Y,'-r v'); 40 ylabel('误差阶数'); 41 title('误差阶数'); 42 toc;
效果图：
二、变步长的隐式Euler方法：
函数文件1：
1 function b=F(t,x0,u,h) 2 b(1,1)=x0(1)-h*x0(2)-u(1); 3 b(2,1)=t*x0(2)+2*h*x0(2)+4*h*t*(2*exp(x0(1))+exp(x0(1)/2))-t*u(2);
函数文件2：
1 function g=Jacobian(x0,t,h) 2 g(1,1)=1; 3 g(1,2)=-h; 4 g(2,1)=4*h*t*(2*exp(x0(1))+0.5*exp(x0(1)/2)); 5 g(2,2)=t+2*h;
函数文件3：
1 function x=Euler(t,u,h) 2 % u：上一节点的数值解或者初值 3 % x0 每次迭代的上一节点的数值 4 % x1 每次的迭代数值 5 % tol 允许误差 6 % f 右端函数 7 x0=u; 8 tol=1e-5; 9 x1=x0-Jacobian(x0,t,h)\F(t,x0,u,h); 10 while (norm(x1-x0,inf)>tol) 11 %数值解的2范数是否在误差范围内 12 x0=x1; 13 x1=x0-Jacobian(x0,t,h)\F(t,x0,u,h); 14 end 15 x=x1;%不动点
脚本文件：
1 tic; 2 clear 3 clc 4 y(1:2,1)=[0;0];%初值 5 e=1e-5;%误差过小 6 tol=1e-5;%指定的误差 7 N=16;%节点的步数 8 h=1/N;%初始步长 9 t=0:h:1; 10 i=1; 11 while t(i)<=1 12 k=1; 13 while k==1 14 y(1:2,i+1)=Euler(t(i)+h,y(1:2,i),h);%符合误差的数值解 15 % y1_half=Euler(h/2,y(1:2,i));%半步长的中点数值解 16 y1_half=Euler(t(i)+h,y(1:2,i),h/2);%半步长的右端点的数值解 17 y1_one=Euler(t(i)+h,y1_half,h/2); 18 Estimate_error=2*norm(y(1:2,i+1)-y1_one);%中间估计误差 19 if Estimate_error<tol%指定误差 20 k=0;%步长相差不大，或者说正好在指定的误差范围内，则确定选择h作为步长。 21 elseif Estimate_error<e%误差过小 22 h=2*h; 23 else%近似估计误差大于指定误差 24 h=h/2; 25 end 26 end 27 t(i+1)=t(i)+h; 28 i=i+1; 29 end 30 f=inline('-2*log(1+x.^2)','x'); 31 Accurate=f(t); 32 subplot(1,3,1) 33 plot(t,y(1,:)); 34 xlabel('t');ylabel('numerical'); 35 title('the image of numerical solution'); 36 grid on ; 37 subplot(1,3,2) 38 plot(t,Accurate); 39 xlabel('t');ylabel('Accurate'); 40 title('the image of Accurate solution'); 41 grid on ; 42 subplot(1,3,3) 43 plot(t,y(1,:)-Accurate); 44 xlabel('t');ylabel('error'); 45 title('the image of error solution'); 46 grid on ; 47 toc;
效果图：
中心差分法：
函数文件1：
1 function b=F(t,x0,h,N) 2 b(1,1)=4*x0(1)-x0(2); 3 b(2,1)=-2*x0(1)+4*h^2*(2*exp(x0(1))+exp(x0(1)/2))+(1+h/t(2))*x0(2); 4 for i=2:N-1 5 b(i+1,1)=(1-h/t(i+1))*x0(i-1)-2*x0(i)+4*h^2*(2*exp(x0(i))+exp(x0(i)/2))+(1+h/t(i+1))*x0(i+1); 6 end
函数文件2：
1 function g=Jacobian(t,x0,h,N) 2 g(1,1)=4; 3 g(1,2)=-1; 4 g(2,1)=-2+4*h^2*(2*exp(x0(1))+1/2*exp(x0(1)/2)); 5 g(2,2)=1+h/t(2); 6 for i=2:N-1 7 g(i+1,i-1)=1-h/t(i+1); 8 g(i+1,i)=-2+4*h^2*(2*exp(x0(i))+1/2*exp(x0(i)/2)); 9 g(i+1,i+1)=1+h/t(i+1); 10 end
函数文件3：
1 function x=newton_Iterative_method(t,u,h,N) 2 % u：上一节点的数值解或者初值 3 % x0 每次迭代的上一节点的数值 4 % x1 每次的迭代数值 5 % tol 允许误差 6 % f 右端函数 7 x0=u; 8 tol=1e-11; 9 x1=x0-Jacobian(t,x0,h,N)\F(t,x0,h,N); 10 while (norm(x1-x0,inf)>tol) 11 %数值解的2范数是否在误差范围内 12 x0=x1; 13 x1=x0-Jacobian(t,x0,h,N)\F(t,x0,h,N); 14 end 15 x=x1;%不动点
脚本文件：
1 tic; 2 clear 3 clc 4 N=[10,20,40,80,160,320,640]; 5 for k=1:length(N) 6 h=1/N(k); 7 x=0:h:1; 8 fun1=inline('-2*log(1+x.^2)'); 9 for i=1:length(x) 10 Accurate(i,1)=fun1(x(i)); 11 end 12 Numerical=zeros(N(k)+1,1); 13 Numerical(2:end,1)=newton_Iterative_method(x,Numerical(2:end,1),h,N(k)); 14 error=Numerical-Accurate; 15 error_inf(k)=norm(error,inf); 16 figure(k) 17 subplot(1,3,1) 18 plot(x,Accurate); 19 xlabel('x'); 20 ylabel('Accurate'); 21 grid on 22 subplot(1,3,2) 23 plot(x,Numerical); 24 xlabel('x'); 25 ylabel('Numerical'); 26 grid on 27 subplot(1,3,3) 28 plot(x,error); 29 xlabel('x'); 30 ylabel('error'); 31 grid on 32 end 33 for i=2:length(N) 34 INF(i-1)=log2(error_inf(i-1)/error_inf(i)); 35 end 36 figure(length(N)+1) 37 plot(1:length(N)-1,INF,'-rh'); 38 xlabel('x'); 39 ylabel('误差阶数'); 40 title('非线性高阶常微分差分格式误差阶'); 41 grid on 42 toc;
效果图：
胡冬冬
