Euler's code for multiple ODEs

7 次查看（过去 30 天）

Nicholas Moung 2021-2-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/739742-euler-s-code-for-multiple-odes

评论： James Tursa 2021-2-8

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

回答（1 个）

James Tursa 2021-2-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/739742-euler-s-code-for-multiple-odes#answer_618267

编辑：James Tursa 2021-2-8

在 MATLAB Online 中打开

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 个评论
显示 10更早的评论隐藏 10更早的评论

Nicholas Moung 2021-2-8

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 2021-2-8

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.

请先登录，再进行评论。

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