solving sets of two differential equations using Runge Kutta 4th order code
85 views (last 30 days)
Show older comments
I would like to know why those two files are not wokring , well , the first one is not running , and the second one is running but giving me wrong solution
I am trying to solve two sets of differential equations using the Runge-kutta 4th order
% Solve y'(t)=-2y(t) with y0=3, 4th order Runge Kutta
y10 = 3; % Initial Condition
y20=1;
h=0.5; % Time step
t = 0:h:2; % t goes from 0 to 2 seconds.
%yexact = 3*exp(-2*t) % Exact solution (in general we won't know this
ystar1 = zeros(size(t)); % Preallocate array (good coding practice)
ysar2=zeros(size(t));
ystar1(1) = y10; % Initial condition gives solution at t=0.
ystar2(1)=y20;
for i=1:(length(t)-1)
k11 = -2*ystar1(i) % Approx for y gives approx for deriv
y11 = ystar1(i)+k11*h/2; % Intermediate value (using k1)
k12 = -2*ystar2(i) % Approx for y gives approx for deriv
y12 = ystar2(i)+k12*h/2; % Intermediate value (using k1)
k21 = -2*y11 % Approx deriv at intermediate value.
y21 = ystar1(i)+k21*h/2; % Intermediate value (using k2)
k22 = -2*y12 % Approx deriv at intermediate value.
y22 = ystar2(i)+k22*h/2; % Intermediate value (using k2)
k31 = -2*y21 % Another approx deriv at intermediate value.
y31 = ystar1(i)+k31*h; % Endpoint value (using k3)
k32 = -2*y22 % Another approx deriv at intermediate value.
y32 = ystar2(i)+k32*h; % Endpoint value (using k3)
k41 = -2*y31 % Approx deriv at endpoint value.
k42 = -2*y32 % Approx deriv at endpoint value.
ystar1(i+1) = ystar1(i) + (k11+2*k21+2*k31+k41)*h/6; % Approx soln
ystar12(i+1) = ystar2(i) + (k12+2*k22+2*k32+k42)*h/6; % Approx soln
end
%plot(t,yexact,t,ystar1);
%legend('Exact','Approximate');
the above code was taken from a site and I was trying to add another equation to use the same code for solving two differentia equations, the error that I keep getting is this one :
Index exceeds the number of array elements (1).
Error in RK4 (line 17)
k12 = -2*ystar2(i) % Approx for y gives approx for deriv
the second file is the folloiwng one ;
% solving sets of differential equations using RK4
%constants
k1=0.001;
k2=0.001;
V1=1;
V2=1;
F=100;
Ca0=0.5;
% initial conditions
t(1)=0;
Ca1(1)=0.5;
Ca2(1)=0.5;
% define function handles
fCa1=@(Ca1, t) (F/V1)*(Ca0-Ca1)-k1*Ca1
fCa2=@(Ca1,Ca2, t) (F/V2)*(Ca1-Ca2)-k2*Ca2
% step size
h=0.1;
tfinal=5;
N=ceil (tfinal/h)
% update loop
for i=1:N
t(i+1)=t(i)+h;
k11=h*fCa1(t(i), Ca1(i));
k12=h*fCa2(t(i), Ca1(i), Ca2 (i));
k21=h*fCa1( t(i) +h/2, Ca1(i)+k11/2);
k22=h*fCa2( t(i) +h/2, Ca1(i)+k11/2, Ca2(i)+k12/2);
k31=h*fCa1( t(i) +h/2, Ca1(i)+k21/2);
k32=h*fCa2( t(i) +h/2, Ca1(i)+k21/2, Ca2(i)+k22/2);
k41=h*fCa1( t(i) +h, Ca1(i)+k31);
k42=h*fCa2( t(i) +h, Ca1(i)+k31, Ca2(i)+k32);
Ca1(i+1)=Ca1(i)+h/6*(k11+2*k21+2*k31+k41);
Ca2(i+1)=Ca2(i)+h/6*(k12+2*k22+2*k32+k42);
end
plot(t,Ca1)
this one is runing put the solution is not correct
I appreicate your help
2 Comments
Answers (1)
Alan Stevens
on 16 Mar 2021
Here is the corrected version 1:
% Solve y'(t)=-2y(t) with y0=3, 4th order Runge Kutta
y10 = 3; % Initial Condition
y20=1;
h=0.5; % Time step
t = 0:h:2; % t goes from 0 to 2 seconds.
yexact = 3*exp(-2*t); % Exact solution (in general we won't know this
ystar1 = zeros(size(t)); % Preallocate array (good coding practice)
ystar2=zeros(size(t));
ystar1(1) = y10; % Initial condition gives solution at t=0.
ystar2(1)=y20;
for i=1:(length(t)-1)
k11 = -2*ystar1(i); % Approx for y gives approx for deriv
y11 = ystar1(i)+k11*h/2; % Intermediate value (using k1)
k12 = -2*ystar2(i); % Approx for y gives approx for deriv
y12 = ystar2(i)+k12*h/2; % Intermediate value (using k1)
k21 = -2*y11; % Approx deriv at intermediate value.
y21 = ystar1(i)+k21*h/2; % Intermediate value (using k2)
k22 = -2*y12; % Approx deriv at intermediate value.
y22 = ystar2(i)+k22*h/2; % Intermediate value (using k2)
k31 = -2*y21; % Another approx deriv at intermediate value.
y31 = ystar1(i)+k31*h; % Endpoint value (using k3)
k32 = -2*y22; % Another approx deriv at intermediate value.
y32 = ystar2(i)+k32*h; % Endpoint value (using k3)
k41 = -2*y31; % Approx deriv at endpoint value.
k42 = -2*y32 ; % Approx deriv at endpoint value.
ystar1(i+1) = ystar1(i) + (k11+2*k21+2*k31+k41)*h/6; % Approx soln
ystar2(i+1) = ystar2(i) + (k12+2*k22+2*k32+k42)*h/6; % Approx soln
end
plot(t,yexact,t,ystar1),grid
legend('Exact','Approximate');
and here is version 2: note my added comments
% solving sets of differential equations using RK4
%constants
k1=0.001;
k2=0.001;
V1=1;
V2=1;
F=100;
Ca0=0.5;
% initial conditions
t(1)=0;
Ca1(1)=1; %0.5;
Ca2(1)=1; %0.5;
% With Ca1(1)= Ca2(1) = 0.5, both fCa1 and fCa2 are zero
% This means neither Ca1 nor Ca2 will change with time
% so all you get are the build-up of numerical errors!
% define function handles
% Make sure the arguments are in the order you are going to call them!
fCa1=@(t, Ca1) (F/V1)*(Ca0-Ca1)-k1*Ca1;
fCa2=@(t, Ca1,Ca2) (F/V2)*(Ca1-Ca2)-k2*Ca2;
% step size
h=0.01;
tfinal=1;
N=ceil(tfinal/h);
% update loop
for i=1:N
t(i+1)=t(i)+h;
k11=h*fCa1(t(i), Ca1(i));
k12=h*fCa2(t(i), Ca1(i), Ca2 (i));
k21=h*fCa1( t(i)+h/2, Ca1(i)+k11/2);
k22=h*fCa2( t(i)+h/2, Ca1(i)+k11/2, Ca2(i)+k12/2);
k31=h*fCa1( t(i)+h/2, Ca1(i)+k21/2);
k32=h*fCa2( t(i)+h/2, Ca1(i)+k21/2, Ca2(i)+k22/2);
k41=h*fCa1( t(i)+h, Ca1(i)+k31);
k42=h*fCa2( t(i)+h, Ca1(i)+k31, Ca2(i)+k32);
Ca1(i+1)=Ca1(i)+1/6*(k11+2*k21+2*k31+k41);
Ca2(i+1)=Ca2(i)+1/6*(k12+2*k22+2*k32+k42);
end
plot(t,Ca1,t,Ca2) ,grid
8 Comments
See Also
Categories
Find more on Boundary Value Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!