Thread Subject:

Simple first order differential equation without dsolve

Hello,

I need to solve this simple differential equation using the Euler method:
Dy=xy
Then I have to plot the approximation and the exact answer on the same graph.

Here is my code so far:


a = 0;
b = 1;
N = input(' How many steps? ');

x = linspace(a, b, N+1);
dx = (b-a)/N;
y(1) = 1;

for i = 1:N
    y(i+1) = x(i)*dx + y(i);
end

figure(1)

plot(x,y,'ok')
hold
fplot('exp(x.^2/2)', [0,1])
title(['N = ', num2str(N)])
xlabel('x')
ylabel('y')




I'm sure ther is something wrong with it because increasing the steps will not increase the accuracy of my approximation.
Using dsolve, I checked the solution which is indeed y=exp(x^2/2)

Please tell me what I am doing wrong.


Jules

On 3/11/2011 11:53 AM, Jules Kassmann wrote:
> Hello,
>
> I need to solve this simple differential equation using the Euler method:
> Dy=xy
> Then I have to plot the approximation and the exact answer on the same graph.
>
> Here is my code so far:
>
>
> a = 0;
> b = 1;
> N = input(' How many steps? ');
>
> x = linspace(a, b, N+1);
> dx = (b-a)/N;
> y(1) = 1;
>
> for i = 1:N
> y(i+1) = x(i)*dx + y(i);
> end
>
> figure(1)
>
> plot(x,y,'ok')
> hold
> fplot('exp(x.^2/2)', [0,1])
> title(['N = ', num2str(N)])
> xlabel('x')
> ylabel('y')
>
>
>
>
> I'm sure ther is something wrong with it because increasing the steps will not increase the
> accuracy of my approximation.
> Using dsolve, I checked the solution which is indeed y=exp(x^2/2)
>
> Please tell me what I am doing wrong.
>
>
> Jules


----------------------------
% solve dy/dx = x*exp(x^2/2) using forward Euler

clear all; close all;

% DATA
a = 0;
b = 1;
N = input(' How many steps? ');
x = linspace(a, b, N+1);
dx = (b-a)/N;

RHS = @(x) x.*exp(x.^2/2);

y = zeros(1,length(x));
y(1) = 1; %initial condition

% algorithm
for i = 1:N
    y(i+1) = y(i) + dx * RHS(x(i));
end

% plot solution
figure(1)

plot(x,y,'ro-');
hold on ;
fplot(@(x) exp(x.^2/2), [x(1),x(end)])
xlabel('x'); ylabel('y');
title(sprintf('solution of dy/dx=x*exp(x^2/2), forward Euler, N=%d',N));
legend('numerical','exact');
---------------------------------


--Nasser

Thank you for your answer Nasser.

What is RHS? Aren't you already using the solution of the differential equation to solve it using the Euler method?



Jules

On 3/12/2011 3:40 AM, Jules Kassmann wrote:
> Thank you for your answer Nasser.
>
> What is RHS? Aren't you already using the solution of the differential equation to solve
> it using the Euler method?
>

RHS is the right hand side of the ode.

when you solve an ODE, such as dy(x)/dx = f(x), then f(x) is given, and I called
this function RHS.

In forward Euler, you write

   (y(i+1)-y(i))/dx = f(x(i))

and then solve for y(i+1)

f(x) is not the solution, y(x) is the solution.

In your original question, you did not tell what the ODE was, you
just plotted the solution at the end.

HTH,
--Nasser

  

"Nasser M. Abbasi" <nma@12000.org> wrote in message <ilfv4i$7l4$1@speranza.aioe.org>...
>
> On 3/12/2011 3:40 AM, Jules Kassmann wrote:
> > Thank you for your answer Nasser.
> >
> > What is RHS? Aren't you already using the solution of the differential equation to solve
> > it using the Euler method?
> >
>
> RHS is the right hand side of the ode.
>
> when you solve an ODE, such as dy(x)/dx = f(x), then f(x) is given, and I called
> this function RHS.
>
> In forward Euler, you write
>
> (y(i+1)-y(i))/dx = f(x(i))
>
> and then solve for y(i+1)
>
> f(x) is not the solution, y(x) is the solution.
>
> In your original question, you did not tell what the ODE was, you
> just plotted the solution at the end.
>
> HTH,
> --Nasser
>
>



But I did tell what the ODE was: Dy=x*y or dy(x)/dx=x*y(x)
The solution for that is y=exp(x^2/2).
I'm sorry but I still don't get it. At first, you don't know that the right hand side is equal to x*exp(x^2/2).

On 3/12/2011 7:05 AM, Jules Kassmann wrote:

>
>
> But I did tell what the ODE was: Dy=x*y or dy(x)/dx=x*y(x)
> The solution for that is y=exp(x^2/2).
> I'm sorry but I still don't get it. At first, you don't know that
> the right hand side is equal to x*exp(x^2/2).

You wrote originally

"fplot('exp(x.^2/2)', [0,1])"

And was comparing the solution you obtained to the above, so I took
that as the solution.

Any way, if your ode is dy/dx = x*y, hen simply make f(x,y) = x*y,
this is the RHS function below.

Mr Euler does not care too much what f(x) is ;)

here it is again, using f(x,y)=x*y

-----------------------------------------
% solve dy/dx = x*y using forward Euler

clear all; close all;

% DATA
a = 0;
b = 1;
N = input(' How many steps? ');
x = linspace(a, b, N+1);
dx = (b-a)/N;

RHS = @(x,y) x*y;

y = zeros(1,length(x));
y(1) = 1; %initial condition

% algorithm
for i = 1:N
     y(i+1) = y(i) + dx * RHS(x(i),y(i));
end

% plot solution
figure(1)

plot(x,y,'ro-');
hold on ;
fplot(@(x) exp(x.^2/2), [x(1),x(end)])
xlabel('x'); ylabel('y');
title(sprintf('solution of dy/dx=x*y, forward Euler, N=%d',N));
legend('numerical','exact');
-------------------------------

--Nasser

Ok, now I understand better how it works now =)
One last thing: now the RHS seems like a function of two variables (x and y) but the whole thing is behaving like it is a single variable function...

Actually, I fully understand now (yes, I do!). It just looks like a function of two variable, but it's not: it's to use in the algorithm for the multiplication of x(i) by y(i).

Thanks a lot for your help Nasser!