CODEWITHSUNDEEP
matlab
user3293979
user3293979

Reputation: 1

Forward Euler Method to solve first order ODEs in Matlab

I wrote this Matlab program that is supposed to solve the IVP du/dx= -5000(u(t) - cos(t)) - sin(t) with u(0)=1. My exact solution should be u(t) = cos(t) but the solution I am getting from Euler's forward in my code is huge in comparison to what it should be and what I calculated but I'm not sure where I've gone wrong in my code. Can you find my error?

function main

dt=5;
u0 = 1;
n=50;

[T,U] = euler(dt, u0, n);

uexact = cos(T);

plot(T,U)

hold on 

plot(T, uexact, 'r')

end

function [T,U]= euler( dt, u0, n)

R= dt/n;
T=zeros(1,n+1);
U=zeros(1,n+1);
U(1)=u0;
T(1) = 0;


for j=1:n   

 U(j+1)= U(j)+ R*(rhs(T(j),U(j)));
 T(j+1)= T(j) + R;
end

end


function dP = rhs(t, P)
P = zeros(1,1);

dP = (-5000)*(P - cos(t)) - sin(t);

end

Upvotes: 0

Views: 749

Answers (1)

sebas
sebas

Reputation: 879

You don not have to set P = zeros(1,1) insde the function that approximate de derivative with the formula.

Moreover, the problem you have is the [numerical unstability of forward Euler method][1]. You need a very small time step to make it converge (because of the large coefficient of P inside function dP).

function main

dt=5;
u0 = 1;
n=50000; % increase the number or subintervals (smaller time step) to get stability

[T,U] = euler(dt, u0, n);

uexact = cos(T);

plot(T,U,'bs')

hold on 

plot(T, uexact, 'r')

end

function [T,U]= euler( dt, u0, n)

R= dt/n;
T=zeros(1,n+1);
U=zeros(1,n+1);
U(1)=u0;
T(1) = 0;

for j=1:n
 % Implicit method converges
 % U(j+1) = ( U(j) - R*(-5000)*cos(T(j)) - R*sin(T(j)))/(1 - R*(-5000));
 U(j+1)= U(j)+ R*(rhs(T(j),U(j)));
 T(j+1)= T(j) + R;
end

end


function dP = rhs(t, P)
%P = zeros(1,1); %% It must not be here

dP = (-5000)*(P - cos(t)) - sin(t);

end


  [1]: http://en.wikipedia.org/wiki/Euler_method

Upvotes: 1

Related Questions