confusion on a new schema for solving ode
조회 수: 1 (최근 30일)
이전 댓글 표시
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
채택된 답변
Bobby Fischer
2021년 1월 15일
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
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!