Results 1 to 6 of 6

Math Help - MatLab IVP Problem

  1. #1
    Newbie
    Joined
    Aug 2009
    Posts
    20

    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
    Attached Files Attached Files
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by exphate View Post
    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
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Aug 2009
    Posts
    20
    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))
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Newbie
    Joined
    Aug 2009
    Posts
    20
    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')
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by exphate View Post
    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
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Newbie
    Joined
    Aug 2009
    Posts
    20
    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
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Matlab LSQ problem
    Posted in the Math Software Forum
    Replies: 2
    Last Post: June 22nd 2010, 02:52 AM
  2. Funny problem. What is cos(pi/2) in matlab?
    Posted in the Math Software Forum
    Replies: 8
    Last Post: January 29th 2010, 10:07 AM
  3. matlab problem
    Posted in the Math Software Forum
    Replies: 7
    Last Post: October 28th 2009, 08:09 AM
  4. Matlab Problem For Help!
    Posted in the Math Software Forum
    Replies: 0
    Last Post: June 16th 2009, 05:15 PM
  5. MATLAB Problem
    Posted in the Math Software Forum
    Replies: 2
    Last Post: January 20th 2009, 06:42 AM

Search Tags


/mathhelpforum @mathhelpforum