Matlab:高阶常微分三种边界条件的特殊解法（中心差分法，高精度导数边界处理）

1 function b=F(f,x0,h,N)
2 %       b(1,1)=x0(1)-h*x0(2)-u(1);
3 %       b(2,1)=x0(2)+h*x0(1)^2-u(2)-h*f;
4 b=zeros(N,1);
5 b(1,1)=4*x0(1)-x0(2);
6 b(2,1)=h^2*x0(1)^2-2*x0(1)+x0(2)-h^2*f(1)
7 for i=3:N
8     b(i,1)=x0(i-2)+h^2*x0(i-1)^2-2*x0(i-1)+x0(i)-h^2*f(i-1);
9 end

 1 function g=Jacobian(x0,h,N)
 2 %   g(1,1)=1;
 3 %   g(1,2)=-h;
 4 %   g(2,1)=2*h*x0(1);
 5 %   g(2,2)=1; 
 6 g=zeros(N);
 7 g(1,1)=4;
 8 g(1,2)=-1;
 9 g(2,1)=2*h^2*x0(1)-2;
10 g(2,2)=1;
11 for i=3:N
12 g(i,i-2)=1;
13 g(i,i-1)=2*h^2*x0(i-1)-2;
14 g(i,i)=1;
15 end

 1 function x=newton_Iterative_method(f,u,h,N)
 2  % u：上一节点的数值解或者初值
 3  % x0 每次迭代的上一节点的数值
 4  % x1 每次的迭代数值
 5  % tol 允许误差
 6  % f 右端函数
 7 x0=u;
 8 tol=1e-11;
 9 x1=x0-Jacobian(x0,h,N)\F(f,x0,h,N);
10 while (norm(x1-x0,inf)>tol)
11     %数值解的2范数是否在误差范围内
12     x0=x1;
13     x1=x0-Jacobian(x0,h,N)\F(f,x0,h,N);
14 end
15 x=x1;%不动点

 1 tic;
 2 clear
 3 clc
 4 N=100;
 5 h=1/N;
 6 x=0:h:1;
 7 y=inline('x.^2.*sin(x).^2+2*cos(x)-x.*sin(x)');
 8 fun1=inline('x.*sin(x)');
 9 Accurate=fun1(x);
10 f=y(x(2:N));
11 Numerical=zeros(N+1,1);
12 Numerical(2:end,1)=newton_Iterative_method(f,Numerical(2:end,1),h,N);
13 error=Numerical-Accurate';
14 subplot(1,3,1)
15 plot(x,Accurate);
16 xlabel('x');
17 ylabel('Accurate');
18 grid on
19 subplot(1,3,2)
20 plot(x,Numerical);
21 xlabel('x');
22 ylabel('Numerical');
23 grid on
24 subplot(1,3,3)
25 plot(x,error);
26 xlabel('x');
27 ylabel('error');
28 grid on
29 toc;