Euler's code for multiple ODEs

44 vues (au cours des 30 derniers jours)
Nicholas Moung
Nicholas Moung le 8 Fév 2021
Commenté : James Tursa le 8 Fév 2021
Hello every one!
May I ask you a favour! I've written a small program for Euler's method to solve multiple ODEs. But I'm encountering an error. Could you please give me a pice of advice. In the following I would like to calculate (x and y). I'm sure I made some mistakes.
Function file for ODEs is
function dydx= FncTG(x,y,Wx_r,Wy_r,Wz_r)
dydx(1)=(Wy_r*sin(y))+(Wz_r*cos(y));
dydx(2)=(Wx_r-(tan (x)*Wy_r*cos(y))-(Wz_r*sin(y)));
dydx= dydx(:);
end
------------------------------------------------------------------------
Euler's method code is
function E = fncEg (f,tspan,y0,x,y,Wx_r,Wy_r,Wz_r)
h = 1/32;
t = 1:h:48;
N = length(t); % A vector to store the time values .
y = zeros (1,N); % Initialize the Y vector .
f=@(x,y)(FncTG);
for i = 1:(N-1)
y (i+1)= y(i)+ h*f(x(i),y(i)); % Update approximation y at t+h
end
E=[t' y'];
end
--------------------------------------------------------------------------
The calculation code is
%Initial values of angles
Theda(1)=0.0025/57.2958;
Gamma(1)=0.0013/57.2958;
h=1/32;
t = 1:h:48;
tspan=[1 48];
y0=[Theda(1),Gamma(1)];
%Importing data from STAT
STAT=importdata('STAT.txt');
nx= STAT(:,10);
ny= STAT(:,11);
nz= STAT(:,12);
Wx_d= STAT(:,7);
Wx_r=deg2rad(Wx_d);
Wy_d= STAT(:,8);
Wy_r=deg2rad(Wy_d);
Wz_d= STAT(:,9);
Wz_r=deg2rad(Wz_d);
%Calculating by Euler method
[x,y]=fncEg(@(x,y)f,tspan,y0,x,y,Wx_r,Wy_r,Wz_r)
I'm looking forward to hearing your adivce. Thank you in advance!

Réponses (1)

James Tursa
James Tursa le 8 Fév 2021
Modifié(e) : James Tursa le 8 Fév 2021
You have two states, so your result should have two values per step. So this code:
y = zeros (1,N); % Initialize the Y vector .
f=@(x,y)(FncTG);
for i = 1:(N-1)
y (i+1)= y(i)+ h*f(x(i),y(i));
should be this instead:
y = zeros (2,N); % Initialize the Y vector . 2 rows per step, not 1
y(:,1) = y0; % initialize first state
%f=@(x,y)(FncTG); % delete this line
for i = 1:(N-1)
y(:,i+1)= y(:,i) + h*f(x(i),y(:,i)); % work with columns of y, not elements of y
And you need to use the proper function handle. So this
[x,y]=fncEg(@(x,y)f,tspan,y0,x,y,Wx_r,Wy_r,Wz_r)
should be this
[x,y]=fncEg(@(x,y)FncTG(x,y,Wx_r,Wy_r,Wz_r),tspan,y0,x,y,Wx_r,Wy_r,Wz_r)
You might consider having your function use a tspan and h that is passed in instead of hard coding this inside the function.
  12 commentaires
Nicholas Moung
Nicholas Moung le 8 Fév 2021
I've used x instead of y(1) and y instead of y(2). If I write the function file for derivatives like above, then will I need to change everywhere, such as
y = zeros (2,N); % Initialize the Y vector . 2 rows per step, not 1
y(:,1) = y0; % initialize first state
%f=@(x,y)(FncTG); % delete this line
for i = 1:(N-1)
y(:,i+1)= y(:,i) + h*f(x(i),y(:,i)); % work with columns of y, not elements of y
and
[x,y]=fncEg(@(x,y)FncTG(x,y,Wx_r,Wy_r,Wz_r),tspan,y0,x,y,Wx_r,Wy_r,Wz_r)
James Tursa
James Tursa le 8 Fév 2021
So, you are going to have to make a choice. Either you work with two different state variables x and y, or you work with one two-element state variable y where y(1) and y(2) are the states. If you work with x and y, then yes you will need to rewrite your propagation equations including both the x and the y propagation code explicitly. If you use the two-element y state vector, then you need to rewrite your derivative equation. You can't have it both ways simultaneously. So pick one method and make your code consistent.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by