4th order Runge Kutta Method Differential System

4th order Runge Kutta Method Differential System. I don't know what's wrong with my code and I have been trying to figure out what went wrong. I'm trying to solve a system with 3 First-Order Differential Equations.

clear all
close all
clc

COne = 1000;
CTwo = 1000;
R = 50;
L = 0.1;
Vt = 0.026;
Is = 1*10^(-8);
a = 0.025/2500

a = 1.0000e-05

b = abs(27.673314419-(2*tan(1.50)))

b = 0.5295

b = 5*sin(2*pi*50*0.017)

b = -4.0451

xprime_func = @(t,x,y,z) ((1/COne)*y-(1/(COne*R))*x);
yprime_func = @(t,x,y,z) ((1/L)*z-(1/L)*x);
zprime_func = @(t,x,y,z) ((1/CTwo)*Is*(exp((abs(10*sin(2*pi*50*t))-z)/(2*Vt))-1)-(1/CTwo)*y);

% Define time interval and step size
tmax=0.025; steps=2500; h=tmax/steps;
% Initial conditions:
x(1)=0; y(1)=0; z(1)=0; t(1)=0;
% Estimate of derivatives and marching in time.
for i=1:steps
 t(i+1)=i*h;

 K(1)=h*xprime_func(t(i),x(i),y(i),z(i));
 L(1)=h*yprime_func(t(i),x(i),y(i),z(i));
 M(1)=h*zprime_func(t(i),x(i),y(i),z(i));


 K(2)=h*xprime_func(t(i)+h/2,x(i)+1/2*K(1),y(i)+1/2*L(1),z(i)+1/2*M(1));
 L(2)=h*yprime_func(t(i)+h/2,x(i)+1/2*K(1),y(i)+1/2*L(1),z(i)+1/2*M(1));
 M(2)=h*zprime_func(t(i)+h/2,x(i)+1/2*K(1),y(i)+1/2*L(1),z(i)+1/2*M(1));

 K(3)=h*xprime_func(t(i)+h/2,x(i)+1/2*K(2),y(i)+1/2*L(2),z(i)+1/2*M(2));
 L(3)=h*yprime_func(t(i)+h/2,x(i)+1/2*K(2),y(i)+1/2*L(2),z(i)+1/2*M(2));
 M(3)=h*zprime_func(t(i)+h/2,x(i)+1/2*K(2),y(i)+1/2*L(2),z(i)+1/2*M(2));

 K(4)=h*xprime_func(t(i)+h,x(i)+1*K(3),y(i)+1*L(3),z(i)+1*M(3));
 L(4)=h*yprime_func(t(i)+h,x(i)+1*K(3),y(i)+1*L(3),z(i)+1*M(3));
 M(4)=h*zprime_func(t(i)+h,x(i)+1*K(3),y(i)+1*L(3),z(i)+1*M(3));

 x(i+1)=x(i)+1/6*(K(1)+2*K(2)+2*K(3)+K(4));
 y(i+1)=y(i)+1/6*(L(1)+2*L(2)+2*L(3)+L(4));
 z(i+1)=z(i)+1/6*(M(1)+2*M(2)+2*M(3)+M(4));

end

plot(t,x,t,y,t,z);

[Torsten]

The Table 7.7 and figure 7.9 is what should my values for my x,y,z and it is insanely way far off.

Please include a mathematical description of your problem with equations and constants used, Table 7.7 and figure 7.9.

I would advice to solve the problem first with a sophisticated MATLAB integrator like ode45 to get a reference solution.

[Prince Nino]

dont mind

"a = 0.025/2500

b = abs(27.673314419-(2*tan(1.50)))

b = 5*sin(2*pi*50*0.017)"

Answer 1

[Torsten]

Note that COne and CTwo are given in muF - thus they should be prescribed as 1000e-6 in your code, I guess.

COne = 1000e-6;
CTwo = 1000e-6;
R = 50;
L = 0.1;
Vt = 0.026;
Is = 1e-8;

xprime_func = @(t,x,y,z) ((1/COne)*y-(1/(COne*R))*x);
yprime_func = @(t,x,y,z) ((1/L)*z-(1/L)*x);
zprime_func = @(t,x,y,z) ((1/CTwo)*Is*(exp((abs(10*sin(2*pi*50*t))-z)/(2*Vt))-1)-(1/CTwo)*y);

f = @(t,x,y,z)[xprime_func(t,x,y,z);yprime_func(t,x,y,z);zprime_func(t,x,y,z)];
F = @(t,u)f(t,u(1),u(2),u(3));

% Define time interval and step size
tmax=0.025; steps=2500; tspan=linspace(0,tmax,steps);
% Initial conditions:
u0 = [0;0;0];
[T,U] = ode15s(F,tspan,u0);

figure(1)
plot(T,U(:,2))
grid on

figure(2)
plot(T,[U(:,1),U(:,3)])
grid on