How to extract variable from function file, and how to check is that result real?

47 views (last 30 days)
I am using matlab function file to solve system of differential equations, and I have variable R in my function file.
function f=fun_z(z,p)
ri=2;
R=ri-z*(ri-1);
sig=1; beta=1;
f=zeros(4,1);
f(1)=-32*1*beta/(R.^4*p(1));
f(2)=(-(2-sig)*8*f(1)/(sig.*R)-f(1)*p(2))/p(1);
f(3)=(-p(2)*f(2)+(2-sig)*(-8*f(2)/R-8*f(1)/(R.*R*p(1)))/sig-f(1)*p(3))/p(1);
f(4)=(-f(2)*p(3)-f(3)*p(2)+(2-sig)*8*(-f(3)/R-f(2)/p(1)+p(2)*f(1)/(p(1).*p(1)))/(sig.*R.*R) -f(1)*p(4))/p(1);
I am not sure did I properly tell to Matlab what is necessary to do with variable R so I extracted it from function file into workspace with command:
[zv,pv,R]=ode45(@fun_z,[1 0],[1; 0; 0; 0])
where I previous in function file changed the first line in:
function f=fun_z(z,p,R)
and added line of code in function file:
global R;
As a result I need to get R where it is dependent on z, so it means R need to have the same number of elements as z, but that is not the case. I cannot see why, but I only get 6 variables? What would be correct way to tell Matlab exactly what R (R is radial coordinate, z is longitudinal coordinate, R is linearly dependent on z, where z is going from 1 to 0, with some step) is? Or how to extract it properly?

Accepted Answer

Jan
Jan on 12 Sep 2018
[zv,pv,R]=ode45(@fun_z,[1 0],[1; 0; 0; 0])
This is a bold and wrong guess. Why do you assume, that ode45 replies an arbitrary global variable from the function to be integrated as 3rd output?
You can either modify the function to accept vectors, then:
function [f, R] = fun_z(z, p)
ri = 2;
R = ri - z * (ri - 1);
sig = 1;
beta = 1;
f = zeros(4, size(p, 2));
f(1, :) = -32 * beta / (R .^ 4 * p(1, :));
f(2, :) = ((-2+sig) * 8 * f(1, :) / (sig .* R) - f(1, :) *. p(2, :)) / p(1, :);
... etc. I cannot test this currently, please debug this by your own.
Then:
[zv, pv] = ode45(@fun_z, [1 0], [1; 0; 0; 0]);
[~, R] = fun_z(zv.', pv.');
Or in your case simply:
[zv, pv] = ode45(@fun_z, [1 0], [1; 0; 0; 0]);
ri = 2;
R = ri - zv * (ri - 1);

More Answers (0)

Products


Release

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by