Matlab:正则Euler分裂

时间：2022-11-10 17:22:27

函数文件1：

function b=F(x0,h,u,tau)

      b(,)=x0()-u();

      b(,)=x0()-u()+*h*1e8*cos(tau)*x0();

函数文件2：

function g=Jacobian(x0,h,tau)

  g(,)=;

  g(,)=;

  g(,)=;

  g(,)=+*h*1e8*cos(tau);

函数文件3：

  function x=Euler(h,x0,u,tau)

  % x0 表示迭代初值

  % u 表示上一节点数值

 tol=1e-;

 x1=x0-Jacobian(x0,h,tau)\F(x0,h,u,tau);

 while (norm(x1-x0,inf)>tol)

     %数值解的2范数是否在误差范围内

     x0=x1;

     x1=x0-Jacobian(x0,h,tau)\F(x0,h,u,tau);

 end

 x=x1;%不动点

脚本文件：

clear

clc

N=[1e2 1e3 1e4 1e5];

for k=:length(N)

h=(pi/)/N(k);

t=:h:pi/;

y(:,)=[;];%原问题初值

u(:,)=y(:,);%子问题一的初值

for i=:N(k)

    u(:,i+)=u(:,i)+h*[-u(,i)^*u(,i)-u(,i)^;u(,i)+1e8*sin(*t(i+))];

    y(:,i+)=Euler(h,u(:,i),u(:,i+),t(i+));

    u(:,i+)=y(:,i+);

end

for i=:N(k)+

    Accurate(:,i)=[cos(t(i));sin(t(i))];

end

error=abs(Accurate-y);

E_max(k)=max(error(,:));

end

% plot(t,Accurate(,:),'-r*')

% hold on

% plot(t,y(,:),'bp')

% grid on

相关文章