Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab and ODE's

  1. Jun 6, 2010 #1
    I'm looking at the following example (http://www.cyclismo.org/tutorial/matlab/control.html):

    y'=x^2-y^2, y(0)=1, where they use the Euler method to approximate numerical solutions.

    This is the code:

    I don't understand the third and fourth line at the top. What is [tex]y = 0*x[/tex] and how did they get y(1) = 1...?

    Last edited: Jun 6, 2010
  2. jcsd
  3. Jun 6, 2010 #2


    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    You want to make sure that the size of the 'y' vector is the same as the size of the domain vector x so you can plot one against the other one.. So they start with y=0*x which should make y a vector of all zeroes that is the same length as x. Setting y(1)=1 means that the first entry of y will be 1. This is the same thing as saying that y, at the value of x that is the first one recorded in the domain vector, should be 1 (so all they've done is chosen an initial condition)
  4. Jun 19, 2010 #3
    Ok. I want to put [tex]y''(x) + \left(\frac{2}{x}\right) y'(x) + y^n = 0[/tex]

    into the following template:

    Code (Text):
    function[y,t] = euler(fun,y0,t0,T,h)
    % [y,t] = euler(fun,y0,t0,T,h) -
    % This function computes the solution to the IVP y'(t) = fun(y,t),
    % y(t0)= y0 , for a given function "fun(t,y)" using Euler's method.
    % The function is defined either via an m-file fun.m, or using the
    % "inline" command. y0 is the initial value, T is the maximum time, h
    % is the stepsize and t0=initial time.
    % The output is a vector containing the approximate solution y_euler.
    y(1) = y0;
    t(1) = t0;
    for i=1:T/h
    y(i+1) = y(i) + h*feval(fun,t(i),y(i));
    t(i+1) = t(i) + h;
    % End of m-file euler.m
    There are a few difficulties. First, I am working with a second order differential equation. Secondly, I am interested in a continuous range of solutions for [tex]n[/tex]. For certain values of n, I have a transformation method that simplifies my equation by eliminating the y'(x). But, the main problem I'm dealing with is coding the "n" part. Whenever I try, an error message in the command box comes up that says 'matrix must be squared' and I don't know what this means.

    I would appreciate any help.

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook