Euler's method for 2nd order ODE

I am trying to run a code for solving 2nd order ODE but I have no idea where is my mistake in the following code:

function f=f(t0,x0,y0,N)
h=0.01;
N=5;
t0=0;
x0=0;
y0=0;
H=zeros(1,N+1);
T=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
H(1)=h; T(1)=t0; X(1)=x0; Y(1)=y0;

for j=1:N
    T(j+1)=T(j)+H(j);
    X(j+1)=X(j)+H(j)*Y(j);
    Y(j+1)=Y(j)+H(j)*ahat(T(j),X(j),Y(j));                      
    H(j+1)=H(j)*(1-H(j));
    
end

plot [T X]
plot [T Y]

[Jon]

Please add further details. In what sense do you think you have a "mistake" ? Is the answer not what you think it should be or are you getting error messages. If you are getting error messages please copy and paste the entire error message to the post.

Also if I were to try and call your function to reproduce the error, what values are you supplying for the input arguments, t0,x0,y0 and N?

[infinity]

Hello,

There are several things that need to clarify to help you solve the problem:

1. What is the error? Are you solving a specific problem? what is the intial data, boundary condition and exact soltuion that you will compare?

2. Can you tell to everyone descrition of the numerical method? I guess it is Euler explicit method. Is it correct?

3. In your function, inputs are t0, x0, y0 and N, but why do you set values of these inputs in your function? It does not expect to do this in functions.

4. the output of the function is f, but you did not assign any valule for it.

[Geoff Hayes]

oday - the output parameter is named f as is the name of your function. You probably don't need the output parameter so could change your function signature to

function f(t0,x0,y0,N)

And, like already stated by Trung, why have input parameters if they are then set to specific values?

Finally, your calls to plot are incorrect. I think that you really want to do

plot(T,X);
plot(T,Y);

instead. Note how the square brackets have been removed.

[oday hazaimah]

Jon:

The error message I am getting is " plot [T X] ".

My inputs as you see above, or you mean they shouldn't be defined there ?

Trung VO DUY:

There are several things that need to clarify to help you solve the problem:

1. The error message I am getting is " plot [T X] ".

2. Yes it is Euler explicit method.

3. I thought I have to set values of inputs in the function, so you mean we cannot do this ?

4. My ahat function is : -2/(1-t)*y+x.

Geoff Hayes:

The reason i am typing the "function f=f(t0,x0,y0,N)" and the command "plot [ ]" the way ther are is because they worked with previous stuff and I thought they are correct .

[infinity]

Hello Oday,

3. I thought I have to set values of inputs in the function, so you mean we cannot do this ?

Of course, you can do this. However, when we write a function, for example f(x,y), we would like to have different value of f as changing values of x and y. In your case, you given specific value of inputs inside the function, therefore you could not change the value of "H, T, X, Y".

4. My ahat function is : -2/(1-t)*y+x.

My question is the output of the function "f(t0,x0,y0,N)". But, it is seen that you do not need the output of this function. So, I suggest you can remove the line

function f=f(t0,x0,y0,N)

sine it does not effect to your code.

1. The error message I am getting is " plot [T X] ".

Now, I guess that you know why do you get this error (since the syntax is not correct in matlab).

[oday hazaimah]

So I have the final draft of the code as follows:

function f(t0,x0,y0,N)

h=0.01;
N=5;
t0=0;
x0=0;
y0=0;
H=zeros(1,N+1);
T=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
H(1)=h; T(1)=t0; X(1)=x0; Y(1)=y0;

for j=1:N
    T(j+1)=T(j)+H(j);
    X(j+1)=X(j)+H(j)*Y(j);
    Y(j+1)=Y(j)+H(j)*ahat(T(j),X(j),Y(j));                      
    H(j+1)=H(j)*(1-H(j));
    
end

plot(T,X);
plot(T,Y);      

But it doesn't give any results nor errors , it just showed some graphs

[Geoff Hayes]

oday - since you initialize

x0=0;
y0=0;

then all of your future updates to X and Y will result in zero values:

X(j+1)=X(j)+H(j)*Y(j);

will be zero since Y will always be zero

Y(j+1)=Y(j)+H(j)*ahat(T(j),X(j),Y(j));  

because ahat is the product of x and y. Try using non-zero initial points for x0 and y0 or you may want to revisit your equations.

[oday hazaimah]

I don't know what is wrong with this code, I have changed as you suggested but neither results no errors.

Are you able to run it in your labtop and see if you get any resutls ?

To avoid calling "ahat '' you may change the 3rd line of the for loop with:

Y(j+1)=Y(j)+H(j)*(-2/(1-T(j))*Y(j)+sin(X(j)));

[Geoff Hayes]

If I change the x0 and y0 to both 1, then the results are non-zero and something appears in the figure (you may need to add a hold on so that the second call to plot retains the graphics object of the previous call to plot).

[oday hazaimah]

you mean 'hold on' between the two plots commands ?

I got the graph but no results.

[Geoff Hayes]

What does "no results" mean? What are the values for T, X, and Y? Use the debugger to determine this.

[oday hazaimah]

I mean no numbers come out when I run the file name with the initials t0,x0,y0,N .

I just see the graph window.

[Geoff Hayes]

no numbers come out...can you clarify where you are expecting to see numbers? In the console window? There are no fprintf's or any other code in your function to write data to the console so not sure what you are expecting to see....

[oday hazaimah]

Geoff,

Usually after we run the code we see the results come out in the command window so I am expecting to see columns with the error to compare between the valuse and then see the graphes of the functions X,Y,T

[Geoff Hayes]

oday - since the code is not explicitly writing anything to the console and because all lines are terminated with a semi-colon, then you will not see any output. As for the error...which array/variable does that correspond to?

Answer 1

[Vaibhav Tomar]

Can you please specify the values of T, Y and X. Without that it's difficult to say why the code is throwing errors. I'd suggest to look for an example for second order ODE algorithm.

You can try this alternate code:

% x1=y
%x2=dy
%then
%dx1=dy=x2
%dx2=d2y=-w^2*y=-w^2*x1
%save this function as yourfcn
function dx=yourfcn(t,x)
w=1
dx(1)=x(2)
dx(2)=-w^2*x(1)
%Then call the ode45 function
[t,x]=ode45(@yourfcn,[0,pi],[0 0])
%Then deduce y
y=x(:,1) 

[oday hazaimah]

what values? aren't these the initial values that we already used above in the code and then the for loop will continue afterwards ?