Compare results of different step size using Euler's method

x1=zeros(10000,1);
x2=zeros(10000,1);
x3=zeros(10000,1);
Time=zeros(10000,1);
for i=1:10001
    t=(i-1)*0.0005;
    Time(i)=t;
    if t==0
        x1(i)=0;
        x2(i)=0;
        x3(i)=2;
    else
        x1(i)=x1(i-1)+x2(i-1)*0.0005;
        x3(i-1)=2-9*sin(x1(i-1))-x2(i-1);
        x2(i)=x2(i-1)+x3(i-1)*0.0005;
end
Xa=x1*180/pi;
plot(Time,Xa,'red')

I can get the plot for one time interval but I want to get it for different time intervals such as at 0.001, 0.01,0.1 and also compare the results on the same plot. how do i do this? please help??

Answer 1

[ME]

I have made a few tweaks to your code and it will now allow you to alter one parameter value (the time step size) and produce a new approximation that will be added to your previous plot.

clear all
h = 0.00001;
x1=zeros(round(5/h),1);
x2=zeros(round(5/h),1);
x3=zeros(round(5/h),1);
Time=zeros(round(5/h),1);
for i=1:(5/h)+1
    t=(i-1)*h;
    Time(i)=t;
    if t==0
        x1(i)=0;
        x2(i)=0;
        x3(i)=2;
    else
        x1(i)=x1(i-1)+x2(i-1)*h;
        x3(i-1)=2-9*sin(x1(i-1))-x2(i-1);
        x2(i)=x2(i-1)+x3(i-1)*h;
    end
end
Xa=x1*180/pi;
hold on
plot(Time,Xa)

Here h is your time step and I have just modified everything else to produce the correct length arrays for storing the results. I also took out the 'red' from your plotting command so that each approximation will automatically plot in a different colour.

[Lucky]

Thanks a lot.. This code works fine... But how do I change the value of h in it.

e.g h = 0.00001;0.005;0.01;

I know how to plot them together but how do I put in different values of h together in one code.

Thanks

[ME]

I guess you'd have to adjust this to something like:

h = sort([0.00001 0.005 0.01],'descend');
for step = 1:numel(h)
    x1=zeros(round(5/h(step)),1);
    x2=zeros(round(5/h(step)),1);
    x3=zeros(round(5/h(step)),1);
    Time=zeros(round(5/h(step)),1);
    for i=1:(5/h(step))+1
        t=(i-1)*h(step);
        Time(i)=t;
        if t==0
            x1(i)=0;
            x2(i)=0;
            x3(i)=2;
        else
            x1(i)=x1(i-1)+x2(i-1)*h(step);
            x3(i-1)=2-9*sin(x1(i-1))-x2(i-1);
            x2(i)=x2(i-1)+x3(i-1)*h(step);
        end
    end
    Xa=x1*180/pi;
    hold on
    plot(Time,Xa)
end

Here the sort is to make sure that you start with the largest step size, otherwise you'll run into issues when you come to plotting.

Answer 2

[Lucky]

Thanks a lot again, but it still doesn't solve my problem.

I need all the plots on a single figure.

clear all
h = sort([0.00001 0.005 0.01 1],'descend');
for step = 1:numel(h)
    x1=zeros(round(5/h(step)),1);
    x2=zeros(round(5/h(step)),1);
    x3=zeros(round(5/h(step)),1);
    Time=zeros(round(5/h(step)),1);
    for i=1:(5/h(step))+1
        t=(i-1)*h(step);
        Time(i)=t;
        if t==0
            x1(i)=0;
            x2(i)=0;
            x3(i)=2;
        else
            x1(i)=x1(i-1)+x2(i-1)*h(step);
            x3(i-1)=2-9*sin(x1(i-1))-x2(i-1);
            x2(i)=x2(i-1)+x3(i-1)*h(step);
        end
    end
end
    Xa=x1*180/pi;
    figure(1)
    subplot(2,1,1)
    plot(Time,Xa)
    hold on
    plot(Time,Xa)
    hold off

[ME]

The code I provided does give all the results on a single figure. Do you mean that you want each in a separate subplot? If so then you could use:

h = sort([0.00001 0.005 0.01],'descend');
for step = 1:numel(h)
    x1=zeros(round(5/h(step)),1);
    x2=zeros(round(5/h(step)),1);
    x3=zeros(round(5/h(step)),1);
    Time=zeros(round(5/h(step)),1);
    for i=1:(5/h(step))+1
        t=(i-1)*h(step);
        Time(i)=t;
        if t==0
            x1(i)=0;
            x2(i)=0;
            x3(i)=2;
        else
            x1(i)=x1(i-1)+x2(i-1)*h(step);
            x3(i-1)=2-9*sin(x1(i-1))-x2(i-1);
            x2(i)=x2(i-1)+x3(i-1)*h(step);
        end
    end
    Xa=x1*180/pi;
    hold on
    subplot(1,numel(h),step)
    plot(Time,Xa)
end

[Lucky]

This works well... Thanks a lot...

Greatly appreciate your help

[Jan Smid]

Hello, may I ask what is the reasoning behind the "5/h(step)" while predefining the correct length arrays? It's not clear to me what exactly does this part of the code do.