confusion on a new schema for solving ode

Question

汉武沈 2021-1-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/717720-confusion-on-a-new-schema-for-solving-ode

评论：汉武沈 2021-1-15

Hello all, recently i've been working on an algorithm for solving ode, called QSS family, which was found by Prof.Ernesto Kofman. I want to run the algorithm on matlab code, but i'm not sure if the code is right. Can anyone help me to tell if it's right?

The code below is an example to solve the ode:

, with

and

all equal 0 at

.

clc
clear
tic
delta_q=1e-4;
q1=0;q2=0;
x1=0;x2=0;
t=0;delta_t=0;tfinal=20;
A=zeros(0,3);
n=0;
while (t<tfinal)
    Dx1=q2;
    Dx2=1-3*q1-4*q2;
    if (Dx1>0)
        delta_x1=delta_q+q1-x1;
    else
        delta_x1=delta_q-q1+x1;
    end 
   if (Dx2>0)
       delta_x2=delta_q+q2-x2;
   else
       delta_x2=delta_q-q2+x2;
   end
       delta_t1=delta_x1/abs(Dx1);
       delta_t2=delta_x2/abs(Dx2);
   if (delta_t1<delta_t2)
       delta_t=delta_t1;
       t=t+delta_t;
       x1=x1+Dx1*delta_t;
       x2=x2+Dx2*delta_t;
       q1=x1;
   else
       delta_t=delta_t2;
       t=t+delta_t;
       x1=x1+Dx1*delta_t;
       x2=x2+Dx2*delta_t;
       q2=x2;
   end
   n=n+1;
   A(n,1)=t;
   A(n,2)=x1;
   A(n,3)=x2;
%    A(n,4)=q1;
%    A(n,5)=q2;
end
A;
 figure(1)
 plot(A(:,1),A(:,2),A(:,1),A(:,3));
toc
       

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bobby Fischer 2021-1-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/717720-confusion-on-a-new-schema-for-solving-ode#answer_598655

在 MATLAB Online 中打开

Hello,

function euler1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
tic
delta_q=1e-4;
q1=0;q2=0;
x1=0;x2=0;
t=0;delta_t=0;tfinal=20;
A=zeros(0,3);
n=0;
while (t<tfinal)
    Dx1=q2;
    Dx2=1-3*q1-4*q2;
    if (Dx1>0)
        delta_x1=delta_q+q1-x1;
    else
        delta_x1=delta_q-q1+x1;
    end
    if (Dx2>0)
        delta_x2=delta_q+q2-x2;
    else
        delta_x2=delta_q-q2+x2;
    end
    delta_t1=delta_x1/abs(Dx1);
    delta_t2=delta_x2/abs(Dx2);
    if (delta_t1<delta_t2)
        delta_t=delta_t1;
        t=t+delta_t;
        x1=x1+Dx1*delta_t;
        x2=x2+Dx2*delta_t;
        q1=x1;
    else
        delta_t=delta_t2;
        t=t+delta_t;
        x1=x1+Dx1*delta_t;
        x2=x2+Dx2*delta_t;
        q2=x2;
    end
    n=n+1;
    A(n,1)=t;
    A(n,2)=x1;
    A(n,3)=x2;
    %    A(n,4)=q1;
    %    A(n,5)=q2;
end
figure(1)
close(1)
figure(1)
subplot(1,3,1)
hold on
axis([0,20 -0.05 0.35])
plot(A(:,1),A(:,2),'b')
plot(A(:,1),A(:,3),'r');
grid on
title('HS')
toc
whos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=linspace(0,20,501);
x0=[0 0];
whos
[t,x]=ode45(@fun,t,x0);
whos
subplot(1,3,2)
hold on
axis([0,20 -0.05 0.35])
plot(t,x(:,1),'k')
plot(t,x(:,2),'k')
grid on
title('edo45')
subplot(1,3,3)
hold on
axis([0,20 -0.05 0.35])
plot(A(:,1),A(:,2),'b')
plot(A(:,1),A(:,3),'r');
plot(t,x(:,1),'k.')
plot(t,x(:,2),'k.')
grid on
title('HS, edo45')
    function [dxdt]=fun(~,x)
        x1p=x(2);
        x2p=1-3*x(1)-4*x(2);
        dxdt=[x1p; x2p];
    end
end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

汉武沈 2021-1-15

Thank you guy! In fact, my original intention was to let others to help me code the algorithm more effeciently, then i felt the code is okay for now. Thank you again anyway！

One more thing, I have upated another question about ode, hope you can help me!

请先登录，再进行评论。