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

24 次查看（过去 30 天）
I G
I G 2018-9-12
回答： Jan 2018-9-12
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?

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2018-9-12
[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);

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by