Euler's code for multiple ODEs
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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!
0 comentarios
Respuestas (1)
James Tursa
el 8 de Feb. de 2021
Editada: James Tursa
el 8 de Feb. de 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 comentarios
James Tursa
el 8 de Feb. de 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.
Ver también
Categorías
Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!