Runge Kutta Fourth order and Interpolation

9 views (last 30 days)
We are trying to solve a system of coupled differential equations by using RK4 for calculating the values of pressures 'p1' and 'p11'. Using the calculated value of pressure we interpolate the value of Energy density 'ed' in every iteration.
The problem we are facing is that the value of pressure 'p11' is not reducing gradually to reach the surface pressure 'p11(end)'.
Kindly help us rectify the mistake in the code.
clear all
close all
clc
a = importdata('rhoce2.dat') ;
%number density a(:,1) ,energy density a(:,2), parallel pressure a(:,3)
% perpendicular pressure a(:,4)
ed(1) = 5.773955099050900e-06; %intial value for energy density 5899
p11(1) = 4.985073388197410e-10; %intial value for parallel pressure
p1(1) = 5.891730398017590e-10 ; %intial value for perpendicular pressure
r(1) = 1e-6 ; %intial value for perpendicular radius
z(1) = 1e-6 ; %intial value for parallel radius
m1(1)= 4*pi*r(1).^2 *z(1) *ed(1) ; %intial value of perpendicular mass
m11(1) = m1(1) ; %intial value of parallel mass
%p11(end) = 2.5e-19 ; % parallel pressure at surface
%p1(end) = 1.3e-15 ; % perpendicu lar pressure at surface
f1 = @(z,r,ed,p11,m11) -(ed+p11)*z*m11.^3 *((r+1)./(((ed+p11)*(r.^2 +z.^2).^6) -(r/p11))) ./(2*pi*(r.^2 + z.^2)) ; %dp11/dz
f2 = @(r,ed) (4/3)*pi*r.^2 *ed ; %dm11/dz
f3 = @(r,z,ed,p1,m1) -4*pi*(ed+p1)*(r.^2+z.^2).^6 *((ed+p1)-p1*(r-m1)) ./(m1.^3 *r.^2) ; %dp1/dr
f4 = @(r,ed,z) (8/3)*pi*r*ed*z ; %dm1/dr
h = 1e-6 ; % step size for perpendicular radius r
j = 1e-6 ; % step size for parallel radius z
i= zeros(1, 15000)
i=1
%{
for i=1:1000
if p11(i) <=1e-20
break
end
%}
while p11(i)>=1e-20
w1 = h* f1(p11(i),z(i),r(i),m11(i),ed(i)) ; %%%% for z component (first equation) p11
w2 = h* f1(p11(i)+0.5*w1,z(i)+0.5*h,r(i),m11(i)+0.5*w1,ed(i)) ;
w3 = h* f1(p11(i)+0.5*w2,z(i)+0.5*h,r(i),m11(i)+0.5*w2,ed(i)) ;
w4 = h* f1(p11(i)+w3,z(i)+h,r(i),m11(i)+w3,ed(i)) ;
x1 = h* f2(ed(i),r(i)) ; %%%%%%%%% for parallel mass ( 2nd equation) m11
x2 = h* f2(ed(i),r(i)) ;
x3 = h* f2(ed(i),r(i)) ;
x4 = h* f2(ed(i),r(i)) ;
y1 = h* f3(p1(i),z(i),r(i),m1(i),ed(i)) ; %%%% for r component (3rd equation) p1
y2 = h* f3(p1(i)+0.5*y1,z(i),r(i)+0.5*j,m1(i)+0.5*y1,ed(i)) ;
y3 = h* f3(p1(i)+0.5*y2,z(i),r(i)+0.5*j,m1(i)+0.5*y2,ed(i)) ;
y4 = h* f3(p1(i)+y3,z(i),r(i)+j,m1(i)+y3,ed(i)) ;
v1 = h* f4(r(i),ed(i),z(i)) ; %%%%%%%% for perpendicular mass(4th equation) m1
v2 = h* f4(r(i),ed(i),z(i)+0.5*j) ;
v3 = h* f4(r(i),ed(i),z(i)+0.5*j) ;
v4 = h* f4(r(i),ed(i),z(i)+j) ;
p11(i+1) = p11(i) + (1/6)* (w1+2*w2+2*w3+w4);
m11(i+1) = m11(i) + (1/6)* (x1+2*x2+2*x3+x4);
p1(i+1) = p1(i) + (1/6)* (y1+2*y2+2*y3+y4);
m1(i+1) = m1(i) + (1/6)* (v1+2*v2+2*v3+v4);
r(i+1) = r(i) + h ;
z(i+1) = z(i) + j ;
%ed(2)= interp1( a(:,2), a(:,4) , p1(2), 'linear');
%ed(i+1)= interp1(a(:,2), p1(i+1));
% ed(i+1)= griddedInterpolant((a(:,2)), (a(:,4)), linear, p1(i+1) );
%F= griddedInterpolant(unique (a(:,4)) , unique ( a(:,2)) );
%ed(i+1) = F(p1(i+1)) ;
%%%interpolation
en = unique ( a(:, 2) ) ; % energy density
pnp = unique ( a(:, 4) ) ; % perpendicular pressure
%p1(2)= 5.790305785983432e-10;
x11= 0; x22= 0; y11= 0; y22= 0;
for k= 2 : 14601
if (p1(i+1) < pnp(k))
x11= en(k-1) ;
x22= en(k);
y11= pnp(k-1);
y22= pnp(k);
break
end
end
ed(i+1) = x11 + (x22-x11) * ((p1(i+1)-y11)/(y22-y11)) ;
i = i+1
end
  2 Comments
Torsten
Torsten on 17 Feb 2023
You have differential equations in r and differential equations in z. You cannot mix them in a "combined" Runge-Kutta-method as you obviously try to do in your code.

Sign in to comment.

Answers (1)

Alan Stevens
Alan Stevens on 17 Feb 2023
Edited: Alan Stevens on 17 Feb 2023
Shouldn't
f1 = @(z,r,ed,p11,m11) -(ed+p11)*z*m11.^3 *((r+1)./(((ed+p11)*(r.^2 +z.^2).^6) -(r/p11))) ./(2*pi*(r.^2 + z.^2)) ;
be
f1 = @(z,r,ed,p11,m11) -(ed+p11)*z*m11.^3 *((r.^2+r)./(((ed+p11)*(r.^2 +z.^2).^6) -(r/p11))) ./(2*pi*(r.^2 + z.^2)) ;
Can't check anything much as you don't supply file rhoce2.dat.
  1 Comment
Arunkarthiheyan
Arunkarthiheyan on 22 Feb 2023
Edited: Arunkarthiheyan on 23 Feb 2023
I checked the code after changing the 'f1' which you have indicated, but still the error persists.
rhoce2 file has been attached herewith.
Kindly help us rectify the problem. Thanks @Alan Stevens

Sign in to comment.

Categories

Find more on Programming 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!