Euler's Method of Numerical Approximations

[kishtik]

Hello,
For Euler's Method of Numerical Approximations, my book (Boyce&DiPrima)
gives this algorithm:

Step 1: define f(t,y)
Step 2: input initial values t0 and y0
Step 3: input step size h and number of steps n
Step 4: output t0 and y0
Step 5: for j from 1 to n do
Step 6: k1 = f(t,y)
y = y + h*k1
t = t + h
Step 7: output t and y
Step 8: end

I tried this on Matlab. My code was:

t=input('Enter t0:');
y=input('Enter y0:');
h=input('Enter step size h:');
n=input('Enter number of steps n:');
k1=input('Define f(t,y):');
result=zeros(n,2);
result(1)=t;
result(1,2)=y;
for j=[2:n+1]
k=k1;
y=y+k*h;
t=t+h;
result(j,1)=t;
result(j,2)=y;
end
disp(' t y');
disp(result);

But I realize that my program calculates k1 only for the initial values... How can I change it so that k1 is actually calculated in every turn of for loop?

Thanks.

 

[LeBrad]

You need to evaluate k = f(t,y) INSIDE the for loop. This way each time you change t and y, you will get a new value for k.

t=input('Enter t0:');
y=input('Enter y0:');
h=input('Enter step size h:');
n=input('Enter number of steps n:');
k1=input('Define f(t,y):');
result=zeros(n,2);
result(1)=t;
result(1,2)=y;
for j=[2:n+1]
k=k1;

k = f(t,y)

y=y+k*h;
t=t+h;
result(j,1)=t;
result(j,2)=y;
end
disp(' t y');
disp(result);

 

[kishtik]

Of course, but the value of k1 is evaluated and assigned at the line
k1=input('Define f(t,y):');
before the loop and the information about the function which should be used to evaluate k's is lost... How can I change it so that it actually evaluates k in the loop using the input k1 as a formula rather than the first of k's?

(Edit: 100th post btw.)

 

[LeBrad]

Of course, but the value of k1 is evaluated and assigned at the line
k1=input('Define f(t,y):');
before the loop and the information about the function which should be used to evaluate k's is lost... How can I change it so that it actually evaluates k in the loop using the input k1 as a formula rather than the first of k's?

(Edit: 100th post btw.)

Ok, I thought your input k1 was just the initial value, but you're saying
k1=input('Define f(t,y):');
is supposed to take in a symbolic function f(t,y)?

 

[kishtik]

Yes, and that symbolic expression should be used to compute the next k's in every turn of the loop.

 

[LeBrad]

Yes, and that symbolic expression should be used to compute the next k's in every turn of the loop.

Ok, I think you want to use inline functions. If you use
k1=input('Define f(t,y):')
and enter something like 't^3+4-5*y' with the ' ', then after that put
f = inline(k1);
and it should define f(t,y) = t^3+4-5*y. Then replace k=k1 with k=f(t,y) insinde your loop. I think that will do what you want.

t=input('Enter t0:');
y=input('Enter y0:');
h=input('Enter step size h:');
n=input('Enter number of steps n:');
k1=input('Define f(t,y):');
f = inline(k1);
result=zeros(n,2);
result(1)=t;
result(1,2)=y;
for j=[2:n+1]
k=f(t,y);
y=y+k*h;
t=t+h;
result(j,1)=t;
result(j,2)=y;
end
disp(' t y');
disp(result);

 

[kishtik]

Thank you! That was just what I was looking for. But now I have another problem...
I modified the code to:

t=input('Enter t0:');
y=input('Enter y0:');
h=input('Enter step size h:');
n=input('Enter number of steps n:');
k1=input('Define f(t,y) in single quotation marks:');
f=inline(k1,'y','t');
result=zeros(n,2);
result(1)=t;
result(1,2)=y;
for j=[2:n+1]
k=f(t,y);
y=y+k*h;
t=t+h;
result(j,1)=t;
result(j,2)=y;
end
disp(' t y');
disp(result);

And tried the example in the book(p.103) with,
t0 = 0,
y0 =1,
h=0.1,
n=4,
f(t,y)=3 + exp(-t) - (1/2)*y,
and got this out:
t _________________y
0________________1.00000000000000
0.10000000000000 1.33678794411714
0.20000000000000 1.65805675248899
0.30000000000000 1.96710763513263
0.40000000000000 2.26609371515479

but the book gave

y0 = 1
y1 = 1.3500
y2 = 1.6730
y3 = 1.9712
y4 = 2.2467

Nearly %2 error for every y. Isn't this too much for an error from rounding or "how my program is constructed" or "how my computer executes arithmetic steps"? The book says "minor variations in the last decimal place" are acceptable, but %2?

 

[LeBrad]

k1=input('Define f(t,y) in single quotation marks:');
f=inline(k1,'y','t');
.
.
.
k=f(t,y);

It looks like you have defined f(y,t) and are evaluating f(t,y). So switch one of the two and it should work.

 

[kishtik]

It does.