MATLAB: Confusion on a new schema for solving ode

solve ode numerically new schema

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

Best Answer

  • 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