2nd order lineair non-homogenous differential equation with EULER not working properly

Hi,

I'm trying to solve n'' + 2n' + 3n = e^t + cos(5t) with EULER to get the same results as GeoGebra CAS has, but it's not showing the same graphs. Also the Matlab doesn't go beyond t=2 although I set it to stop at tf=5.

Can someone see what's wrong?

Thank you a lot!

Matlab:

% n''+2n'+3n=0, n(0)=0, n'(0)=2

%Defining functions

first=@(n,x,t) x;

second=@(n,x,t) exp(t) + cos(5*t) - 3*n - 2*x;

%step size

T=.000001;

%max t value

tf=5; %doesn't plot beyond tf=2

%Initial conditions

t(1)=-3;

n(1)=0; %n(0)=0

n2(1)=2; %n'(0)=2

%euler approximation

for i=1:(tf/T)

t(i+1)=t(i)+T;

n(i+1)=n(i)+T*first(n(i),n2(i),t(i));

n2(i+1)=n2(i)+T*second(n(i),n2(i),t(i));

end

plot(t,n)

GeoGebra CAS:

Answer 1

[Torsten]

The number of steps to take to reach tend is (tend-tstart)/dT, not tend/dT.

Further, you apply the initial conditions at t = -3, not t = 0.

[Manon]

thank you. How do I apply initial conditions at t=0 but plot graph from -3 to 3?

[John D'Errico]

@Manon

Just solve the problem TWICE. ODE45 is an initial value solver. So you need to specify the initial conditions at one end of the interval (Usually the left end for some reason, but that is not important here.)

You cannot specify the initial condition at a point inside the region, and then solve the problem out in both directions at once. However, nothing stops you from solving it twice. For example, with

dy/dt = sin(t) + 2

with the initial condition of

y(0)=1.

and an initial condition at t == 0.

But suppose we wish to plot the solution from -3 to 3? Yes, we could use dsolve to solve the problem symbolically. Then we could do anything we want. But how to solve it using ODE45?

dyfun = @(t,y) sin(t) + 2;
[tout1,yout1] = ode45(dyfun,[0,3],1);
[tout2,yout2] = ode45(dyfun,[0,-3],1);

Now, plot the two solutions.

plot(tout1,yout1,'r-',tout2,yout2,'b-',0,1,'go')
xline(0,'g')
grid on

As you can see, they join up at t==0.

[Torsten]

% n''+2n'+3n=0, n(0)=0, n'(0)=2
%Defining functions
syms n(x)
Dn = diff(n,x);
eqn = diff(n,x,2)+2*Dn+3*n==exp(x)+cos(5*x);
conds = [Dn(0)==2,n(0)==0];
nsol(x) = dsolve(eqn,conds);
nsol = simplify(nsol)

nsol(x) = 

f = matlabFunction(nsol);
first=@(n,x,t) x;
second=@(n,x,t) exp(t) + cos(5*t) - 3*n - 2*x;
%step size
T=.01;
%max t value
tstart = 0.0;
tend = 3.0; 
%Initial conditions
tf(1)=tstart;
nf(1)=0; %n(0)=0
nf2(1)=2; %n'(0)=2
%euler approximation
for i=1:(tend/T)
    tf(i+1)=tf(i)+T;
    nf(i+1)=nf(i)+T*first(nf(i),nf2(i),tf(i));
    nf2(i+1)=nf2(i)+T*second(nf(i),nf2(i),tf(i));
end
tstart = 0.0;
tend = -3.0;
tb(1) = tstart;
nb(1) = 0.0;
nb2(1) = 2.0;
%euler approximation
for i=1:(-tend/T)
    tb(i+1)=tb(i)-T;
    nb(i+1)=nb(i)-T*first(nb(i),nb2(i),tb(i));
    nb2(i+1)=nb2(i)-T*second(nb(i),nb2(i),tb(i));
end
hold on
plot([tb(end:-1:2),tf(1:end)],[nb(end:-1:2),nf(1:end)])
plot([tb(end:-1:2),tf(1:end)],f([tb(end:-1:2),tf(1:end)]))
hold off
grid on

[Manon]

thank you a lot Torsten! I've been at it all day yesterday and couldn't get anywhere.

It's weird because I worked the problem out on paper as well and now Matlab gives another awnser to the equations as GeoGebra CAS does...

I worked out C1 and C2 with the given conditions and that gives me the awnser below.

Weird... :(

[Torsten]

I also plotted the "solution" from GeoGebra CAS and there are great differences compared to your Euler solution, especially for x < 0. As you can see, MATLAB solution and Euler solution are almost identical.

My guess is that you incorrectly incorporated the two initial conditions for n and n' at x=0.

But it's no problem: if you insert the two "solutions" into the differential equation, you can easily see which of them is correct.