# Thread: MatLab IVP Problem

1. ## MatLab IVP Problem

Hey guys!

I'm working through some book questions for my numerical methods subject, and I'm having trouble with this.

The IVP is
$(1+t^2) dy/dt + 2ty = -pi*sin(t*pi)$

Now, I got a stack of code such that I could plot the comparison of the numerical and analytical solutions out, and my tutor agreed. (I'll attach it in a word document).

But I'm having issues with finding the "(relative) truncation error" at various values for h. Now, I got out some answers that didn't 'seem right' (based on observation from the different graphs) but yeh. I'm wondering if anyone is able to throw some ideas or even code my way to help me get this out?

It'd be MASSIVELY appreciated

2. Originally Posted by exphate
Hey guys!

I'm working through some book questions for my numerical methods subject, and I'm having trouble with this.

The IVP is
$(1+t^2) dy/dt + 2ty = -pi*sin(t*pi)$

Now, I got a stack of code such that I could plot the comparison of the numerical and analytical solutions out, and my tutor agreed. (I'll attach it in a word document).

But I'm having issues with finding the "(relative) truncation error" at various values for h. Now, I got out some answers that didn't 'seem right' (based on observation from the different graphs) but yeh. I'm wondering if anyone is able to throw some ideas or even code my way to help me get this out?

It'd be MASSIVELY appreciated
Lets assume that you want the relative error at each of the time points for each of your time steps.

Let the numerical solution be computed at $t_1,\ ..,\ t_n$ and have values $yN_1,\ ..,\ yN_n$, and the analytic solution at the same time points be $yA_1,\ ..,\ yA_n$

Then the relative error at $t_i$ is:

$\varepsilon_i=\frac{|yN_i-yA_i|}{|yA_i|}$

or in Matlab like code:

Code:
err=abs(yN-yA)./abs(yA);
The key point is the analytic solution is evaluated at the same points as the numerical.

CB

3. Ive got something similar to that now. BUT, I have the following errors:-
h= 0.1: 0.06197684648065
h= 0.01: 0.00059373260547
h= 0.001: 0.00000593482743
h= 0.0001: 0.00000005934802

I'm worried about the first one (h=0.1) because it doesn't match the pattern that the questions are talking about (O(h^2))

4. Code being:-

%Lab Week 4

%setup DE
f = @(t,y) (-pi*sin(pi*t)-(2*t*y)/(1+t^2));

%int. conditions
a = 0;
b = 1;

xEnd = a + 5; % integrate two units from starting value

format long
%Setup analytic solution

yExact = @(t) (-cos(pi*t)/(1+t^2))

%set up some step sizes to test
h1 = 0.1;
h2 = h1 / 2;
h3 = h2 / 2;
h4 = h3 / 2;
h5 = h4 / 2;
h6 = h5 / 2;

%local truncation error requires one step
[f1, t1] = forwardEuler(f, a, b, h1, 2);
[f2, t2] = forwardEuler(f, a, b, h2, 2);
[f3, t3] = forwardEuler(f, a, b, h3, 2);
[f4, t4] = forwardEuler(f, a, b, h4, 2);

err1 = (abs((( f1(2)-yExact (t1(2)))))/abs(yExact (t1(2))))
err2 = (abs((( f2(2)-yExact (t2(2)))))/abs(yExact (t2(2))));
err3 = (abs((( f3(2)-yExact (t3(2)))))/abs(yExact (t3(2))));
err4 = (abs((( f4(2)-yExact (t4(2)))))/abs(yExact (t4(2))));

stepSize = [h1, h2, h3, h4, h5, h6]; %step sizes into array for plotting
LTErrors = [err1, err2, err3, err4, err5, err6]; %put all errors into an array

plot(stepSize, LTErrors, 'b.', 'linewidth', 2)

xlabel('step size')
ylabel('local truncation error')

5. Originally Posted by exphate
Ive got something similar to that now. BUT, I have the following errors:-
h= 0.1: 0.06197684648065
h= 0.01: 0.00059373260547
h= 0.001: 0.00000593482743
h= 0.0001: 0.00000005934802

I'm worried about the first one (h=0.1) because it doesn't match the pattern that the questions are talking about (O(h^2))
That is nothing to worry about, the $O(h^2)$ behaviour is an asymptotic thing. When $h$ is relativity large other effects can can be seen (such as instability and higher order terms in the error)

CB

6. Cheers dude. I spoke to one of the guys at uni today (masters student) and he went through it in person, and pointed out that the 'large' step size often wont follow the same number pattern.

Thanks for the clarification