Euler's code for multiple ODEs

2021 2 8
1 个回答
更新时间：2021 2 8
7 次查看（30 天）

0 个投票

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!

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

回答（1 个）

James Tursa
James Tursa 2021-2-8
编辑：James Tursa 2021-2-8

1 个投票

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
Nicholas Moung 2021-2-8
编辑：Nicholas Moung 2021-2-8
I did like the following.
function E = fncEg (f,tspan,y0,N)
h = 1/32;
t = 1:h:48;
N = length(t); % A vector to store the time values .
y = zeros (2,N); % Initialize the Y vector . 2 rows per step, not 1
y(:,1) = y0; % initialize first state
f = @(x,y)FncTG(x,y,Wx_r,Wy_r,Wz_r); % use the proper function handle
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
end
And I still have error in output like the following. Because I use this function [x,y]=fncEg(@(x,y)ff,tspan,y0) for calculation x and y.
Error using fncEg
Too many output arguments.
Error in Euler (line 25)
[x,y]=fncEg(@(x,y)ff,tspan,y0)
Thank you Sir! I'm counting on you Sir.
James Tursa
James Tursa 2021-2-8
You can't define the function handle f inside the fncEg function because it can't see the Wx_r,Wy_r,Wz_r variables. You are passing in the f as an argument. So delete this line as I have shown.
Also, it looks like you will need to fix up your derivative function to handle a 2-element y vector as an input. Looks like it is not set up for that properly.
What is the differential equation you are trying to solve?
Nicholas Moung
Nicholas Moung 2021-2-8
The differential equations are
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
James Tursa
James Tursa 2021-2-8
No. That is the code you wrote for the differential equation. I am asking you to post the actual differential equation. E.g., an image from your assignment.
Nicholas Moung
Nicholas Moung 2021-2-8
In this case, Theda is x and Gamma is y Sir. And Wx,Wy,Wz arrays.
James Tursa
James Tursa 2021-2-8
What is the tgv part of this? Is t supposed to be time, and g supposed to be gravitational constant? Why do you have tan(x) in your code?
Nicholas Moung
Nicholas Moung 2021-2-8
In the second eq, tgv is tan (v) which is theda. And I've used x instead of v in my code Sir.
James Tursa
James Tursa 2021-2-8
编辑：James Tursa 2021-2-8
This does not look like tan(v) to me. It looks like it should be t*g*v. Why do you think it should be tan(v)? Where does the tangent function come from?
Assuming it really is tan(v) and this is just some funny notation, you will need to alter your derivative function to use y(1) and y(2) instead of x and y.
Nicholas Moung
Nicholas Moung 2021-2-8
In Russia, we write tg instead of tan. And in this case v is not v but theda which is angle of pitch of aircraft Sir. I still get error. I think I did some mistakes in creating the function of ODEs, didn't I?
Undefined function or variable 'x'.
Error in Euler (line 27)
[x,y]=fncEg(@(x,y)ff(x,y,Wx_r,Wy_r,Wz_r),tspan,y0,x,y,Wx_r,Wy_r,Wz_r)
James Tursa
James Tursa 2021-2-8
If it is tan(v), then your derivative function should look something like this based on the image
function dydx= FncTG(x,y,Wx_r,Wy_r,Wz_r)
dydx(1) = Wy_r*sin(y(2)) + Wz_r*cos(y(2));
dydx(2) = Wx_r - tan(y(1))*(Wy_r*cos(y(2))-Wz_r*sin(y(2)));
dydx= dydx(:);
end
Nicholas Moung
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
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by