Euler's code for multiple ODEs

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!

Answer 1

[James Tursa]

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.

[Nicholas Moung]

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]

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]

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]

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]

In this case, Theda is x and Gamma is y Sir. And Wx,Wy,Wz arrays.

[James Tursa]

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]

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]

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]

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]

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]

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]

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.