How to solve an equation with arrays
17 views (last 30 days)
Show older comments
I currently have this code where "eqn" has to be solved to find k_eff and then b vs k_eff must be plotted. This is what I have so far and it seems like everythibng is just iterating in an infinite loop. I'm not sure what's wrong, so any help would be appreciated.
clc
clear all
close all
% reflector
Dr=0.088;
Sa_r=0.0033;
%core
Dc=1.45;
Sa_c=0.0044;
nu=2.5;
Sf=0.0032;
b=0.01:0.1:500;
a=30;
syms k_eff
k_eff=sym('k_eff',length(b));
for i=1:length(b)
for j= 1:length(k_eff)
eqn(j,i)=Dc*sqrt(-(2.405/(b(i)/2+2*Dc))^2+((nu*Sf)/k_eff(j)-Sa_c)/Dc)*tanh(sqrt((2.405/(b(i)/2+2*Dr))^2+Sa_r/Dr)*(a+2*Dr))-Dr*sqrt((2.405/(b(i)/2+2*Dr))^2+Sa_r/Dr)*tan(sqrt(-(2.405/(b(i)/2+2*Dc))^2+((nu*Sf)/k_eff(j)-Sa_c)/Dc)*(-b(i)-2*Dc));
sol=solve(eqn,k_eff(j))
plot(b(i),k_eff(j))
hold
on
grid on
end
end
plot(b(i),k_eff(j))
hold
on
grid on
0 Comments
Accepted Answer
Torsten
on 13 Feb 2023
Edited: Torsten
on 13 Feb 2023
syms b k_eff
% reflector
Dr=0.088;
Sa_r=0.0033;
%core
Dc=1.45;
Sa_c=0.0044;
nu=2.5;
Sf=0.0032;
a=30;
eqn = Dc*sqrt(-(2.405/(b/2+2*Dc))^2+((nu*Sf)/k_eff-Sa_c)/Dc)*tanh(sqrt((2.405/(b/2+2*Dr))^2+Sa_r/Dr)*(a+2*Dr))-Dr*sqrt((2.405/(b/2+2*Dr))^2+Sa_r/Dr)*tan(sqrt(-(2.405/(b/2+2*Dc))^2+((nu*Sf)/k_eff-Sa_c)/Dc)*(-b-2*Dc)) == 0;
B=0.01:1:500;
k_eff0 = 0.01;
for i = 1:numel(B)
K_eff(i) = double(vpasolve(subs(eqn,b,B(i)),k_eff0));
k_eff0 = K_eff(i);
end
plot(B,K_eff)
grid on
More Answers (1)
William Rose
on 13 Feb 2023
The commands
b=0.01:0.1:500;
syms k_eff
k_eff=sym('k_eff',length(b));
cause Matlab to try to create a sybolic array k_eff with dimensions 5000x5000. That alone will take a very long time and a lot of memory. I dount that is really what you want. What DO you want k_eff to look like?
Also, you have a for loop inside another for loop, and the inner loop has a plot command. This would lead to the generation of 25 million plots. No wonder the code seems stuck in an infinite loop. I suggest you figure out what you really want to do. Try it first with b=0.01:.1:1. Is it doing what you want?
Good luck.
0 Comments
See Also
Categories
Find more on Optimization 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!