Evaluating a 2nd order ODE using the Runge-Kutta method

46 Ansichten (letzte 30 Tage)
Greetings,
I've been working on a 2nd order ODE: y''(t) = -e^(3t)*y'(t) - y(t) + (5-2e^(-3t))*e^(-2t) +1
With initial conditions y(0) = 2 ; y'(0) = -2 where 0<t<1
I need to use a third order Runge Kutta method, which I can code for a 1st order ODE.
However given a 2nd order differential equation, I'm having difficulties implementing the ODE into my Runge Kutta code.
I've tried writing the 2nd order ODE in a linear system of 1st order ODEs but im still stuck.
By the way, i do not want to use ode45 or ode23 commands. Thank you for your feedback!
Here is what I've been working on so far:
function RK3 = func(a,b,n)
a = 0;
b = 1;
n = input('Enter the number of intervals:');
h = (b-a)/n;
y(1) = 2;
y(2) = -2;
for i=1:n
t = (i-1)*h;
k1 = f(t,y(:,i));
k2 = f(t+0.5*h,y(:,i)+0.5*k1*h);
k3 = f(t+h,y(:,i)-k1*h+(2*k2*h));
y(i+1)= y(:,i)+h/6*(k1+4*k2+k3); %approximated value of the ODE using RK3
end
for k=1:n+1
t=a+(k-1)*h;
exact(k)=exp(-2*t)+1; %exact value of the ODE
end
error = max(abs(y-exact));
max(exact)
max(y)
error
end
function fty = f(t,y) %matrix form of the 2nd order ODE
y = [0 1; -1 -exp(-3*t)];
u = [0; 5-2*exp(-3*t)*exp(-2*t)+1];
fty = y + u;
end

Akzeptierte Antwort

Alan Stevens
Alan Stevens am 31 Dez. 2020
More like this:
a = 0;
b = 1;
%n = input('Enter the number of intervals:');
n = 100;
h = (b-a)/n;
y = [2; -2];
t = 0:h:n*h;
for i=1:n
k1 = f(t(i),y(:,i));
k2 = f(t(i)+0.5*h,y(:,i)+0.5*k1*h);
k3 = f(t(i)+h,y(:,i)-k1*h+2*k2*h);
y(:,i+1)= y(:,i)+h/6*(k1+4*k2+k3); %approximated value of the ODE using RK3
end
exact = exp(-2*t)+1;
disp(max(abs(exact-y(1,:))))
plot(t,y(1,:),'r',t,exact,'b--'),grid
xlabel('t'),ylabel('y')
legend('RK3','Exact')
function fty = f(t,y) %matrix form of the 2nd order ODE
Y = [0 1; -1 -exp(-3*t)]*y;
u = [0; (5-2*exp(-3*t))*exp(-2*t)+1];
fty = Y + u;
end
  1 Kommentar
Emir Alp Karslioglu
Emir Alp Karslioglu am 1 Jan. 2021
This code works like a charm! Thank you so much for your assistance. I need to work more on matlix manipulation in matlab it seems. Much appreciated!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by