How to solve an equation with arrays

13 views (last 30 days)
Din N
Din N on 13 Feb 2023
Edited: Torsten on 13 Feb 2023
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

Accepted Answer

Torsten
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
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.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!