how to iterate an ODE

I need to run the code written in the first following section for the variable "P" running from 400000 to 900000, then consider only the max( sol.y) value considered in the sol.x>50 and sol.x<180 range for every "P(i)" solution.

i've tried to fit the code in a for loop i=1:length(n) allocating a P(i) value to the P variable where n is a "linspace" vector but it is not correct, it says that is not possible to convert from sym to double.

Could someone help me please?

this is the code whitout the loop:

%variables:

P=0.5*800000;%kN                                -carico agente
R=20;%mm                                        -raggio medio
L=3*75;%mm                                      -lunghezza tubo
s=1;%mm                                         -spessore di parete
ni=0.3;%                                        -coeff poisson
E=210000;%MPa                                   -modulo young
w=2;%MPa                                        -pressione int/ext
K=0.5;%                                         -incasto incastro



E1=E/(1-(ni^2));%                     -modulo young dilataz laterale impedita
I=(2*pi*R*(s^3))/12;%                 -momento di inerzia cilindro
k=2*pi*((E*s)/R);%                    -rigidezza molla fondazione
Pcr=sqrt(4*k*E1*I)%                   -carico critico
Pcre=((pi)^2*E*I)/(K*((L/1000)^2))%   -carico critico secondo Eulero

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%ODE:

f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I))*((P*u(3))+(k*u(1)))]; 

bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001]; 


zmesh = linspace(0, L, 100000);

iguess= @(z) [0; 0; 0; 0];

solinit = bvpinit(zmesh,iguess ); 
opts = bvpset('RelTol',0.1,'AbsTol',0.1,'Stats','on'); 
                                                       
sol= bvp4c(f, bc, solinit,opts);


hold on
grid on
plot(sol.x(1,:),sol.y(1,:))


legend('w=0')
xlabel('z');
ylabel('solution u');

this is the for loop that should solve the differential equation f for every P value:

syms P;
n=linspace(400000,900000);

for i=1:length(n)
f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I))*((P(i)*u(3))+(k*u(1)))];
bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
zmesh = linspace(0, L, 100000);
iguess= @(z) [0; 0; 0; 0];
solinit = bvpinit(zmesh, iguess);
Sol= bvp4c(f, bc, solinit);
end

the section that finds the max of sol.y is:

x = sol.x;
y = sol.y;
idx = (x > 50) & (x < 180);
x_new = x(idx);
y_new = y(1, idx);

y_new_max = max(y_new)

Answer 1

[Walter Roberson]

i've tried to fit the code in a for loop i=1:length(n) allocating a P(i) value to the P variable where n is a "linspace" vector but it is not correct, it says that is not possible to convert from sym to double.

Let us have a look at the code:

syms P;

That code is equivalent to

P = sym('P')

which stores a reference to the symbolic variable P (that lives inside the symbolic engine) in the MATLAB variable P . P at the MATLAB level is a scalar -- a reference to exactly one symbolic variable.

f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I))*((P(i)*u(3))+(k*u(1)))];

When an anonymous function is built, nothing in the function body is executed. The exact code is copied, using the following rules:

1. Any variable listed in the @() block (that is not shadowed by further @ levels) is detected inside the block of code, and replaced by a reference to the corresponding positional parameter number. So for example the u(2) part would be replaced by code that says "at execution time, index the second positional parameter passed in, at linear index 2"
2. Any variable that is not listed in the @() is searched for in the current local environment -- local variable assignments already executed, parameter names in the function definition, assignments already executed in any function that the current function is nested in. If found, then the current value of that variable is copied in to the data structure for the anonymous function, and the variable name will be replaced by a reference to that copy of the data. After the function is built, no changes to that variable outside of the anonymous function execution can change the value of the variable for the purposes of the anonymous function.
3. Any variable that is not listed in the @() and not found in the above search is marked as undefined. In current releases of MATLAB, when the anonymous function is executed and there is a variable marked as undefined, there will be an error message generated at the point the value of the variable is needed. This is a change from some of the older versions of MATLAB: in some of the older versions of MATLAB, if any variable was undefined after the search, then all variables except the parameters were marked as needing to be resolved, potentially leading to a dynamic run-time resolution of values instead of a captured-value evaluation.

Notice in case 2 that any indexing is not done at the time the function handle is built: the entire value of the variable is captured and the indexing is done at run-time.

Putting these together: your P is scalar (a reference to a single symbolic variable). When you build the anonymous function f, no matter what the value of the index i is, the entire (scalar) P is copied into the function handle f. Then when you execute the function handle, the scalar symbolic variable P will be indexed by the value of i that was captured.

On the first iteration, when i is 1, then it is valid to index the scalar symbolic variable P at 1, giving back the scalar symbolic variable. That symbol value then gets used and since it is symbolic, that triggers all of the other items in the [] to be converted to symbolic, giving you a symbolic output from invoking f. But the boundary value functions and ode*() require that the returned value be strictly single or double. Even if somehow P were assigned a symbolic number so that the output could be converted from symbolic to integer, you would later have problems with ode45 complaining that the result was symbolic (symbolic numbers are still symbolic datatype.) But before you get there, the bvp4c routine tries to assign the result to a numeric array, and since you have the unresolved symbolic variable P in the expression, that is not going to work.

If you somehow got through that to the second iteration when i is 2, then you would be indexing the scalar symbolic variable P at location 2, which would fail.

All of which is to say that you are not allocating a P(i) value to the P variable like you claim you are. You are not assigning to P(i) anywhere; the only assignment to P is the

syms P

Perhaps instead of using syms P you should be assigning P as a numeric vector, possibly the result of doing a linspace() or a calculation processing a linspace() result.

[federico midei]

thank u !