Using a function and fsolve with a for loop

2 views (last 30 days)
I'd like to solve the system of 4 non linear equations and 4 unknonws over a span of values for T.
I'm recieving the following error:
"This statement is not inside any function.
(It follows the END that terminates the definition of the function "eqn".)"
function F = eqn(x,T)
F(1) = 10.*(T-x(1))+0.2.*(T.^4-x(1).^4) - x(4);
F(2) = -1.*(1/x(1)-1/x(2)) - x(4);
F(3) = 2.*(x(3)-x(2))-x(4);
F(4) = 5.*(1-x(3))-x(4);
end
T = 1:1:40;
for i = 1:numel(T)
x(i,:)=fsolve(@(x)eqn(x,T(i)),[0,0,0,0]);
end

Answers (2)

Steven Lord
Steven Lord on 10 Feb 2023
Either remove the section of code starting with "T = 1:1:40;" from your function file and run it in the Command Window or another file, or move that section of code to before the line "function F = eqn(x,T)". This second approach may require you to change the name of your function or the name of the file that contains it.
  3 Comments
Torsten
Torsten on 10 Feb 2023
Yes, you divide by 0 in the expressions 1/x(1) and 1/x(2) if you set [0 0 0 0] as initial guess.
Alex Sha
Alex Sha on 11 Feb 2023
See the results below:
T x1 x2 x3 x4
1 1 1 1 -4.41868100791151E-15
2 -0.04221805467944030 -15.5355259373919 -3.72443598211183 23.6221799105584
3 -0.02152989066740660 -31.4907092033764 -8.28305977239318 46.4152988619656
4 -0.01094985711415290 -62.9166489977652 -17.2618997136447 91.309498568229
5 -0.00571215265706289 -121.539985068469 -34.0114243052774 175.057121526389
6 -0.00313248055977908 -222.461927363911 -62.8462649611353 319.231324805659
7 -0.00181745226754944 -384.152722165812 -109.043634904507 550.218174522591
8 -0.00111208392250994 -628.447784587467 -178.842224167857 899.211120839266
9 -0.000713160880548333 -980.544992126182 -279.441426321773 1402.2071316088
10 -0.000476189242038096 -1469.00333332469 -419.000952378488 2100.00476189245
11 -0.000329141847745449 -2125.74230399299 -606.640658283718 3038.20329141854
12 -0.000234345559697521 -2986.04164041889 -852.440468691099 4267.20234345554
13 -0.000171168338151485 -4088.54119817842 -1167.44034233668 5842.20171168344
14 -0.000127824907104757 -5475.2408947744 -1563.64025564984 7823.20127824912
15 -9.73235904377666E-5 -7191.50068126512 -2054.00019464718 10275.0009732358
16 -7.53738494234600E-5 -9286.04052761695 -2652.44015074771 13267.2007537385
17 -5.92620663416646E-5 -11810.9404148345 -3373.84011852413 16874.2005926207
18 -4.72250545218803E-5 -14821.6403305754 -4234.04009445011 21175.2004722506
19 -3.80891432006883E-5 -18376.940266624 -5249.84007617831 26254.2003808915
20 -3.10559002788023E-5 -22539.0002173913 -6439.00006211179 32200.000310559
21 -2.55713925803802E-5 -27373.3401789998 -7820.2400511428 39106.200255714
22 -2.12444126098758E-5 -32948.8401487108 -9413.24004248879 47071.2002124441
23 -1.77941641615093E-5 -39337.7401245591 -11238.6400355883 56198.2001779416
24 -1.50160972175629E-5 -46615.6401051127 -13318.0400300322 66595.2001501609
25 -1.27591706301686E-5 -54861.5000893142 -15674.0000255183 78375.0001275917
26 -1.09104556930710E-5 -64157.6400763732 -18330.0400218209 91655.2001091045
27 -9.38454289762297E-6 -74589.7400656918 -21310.6400187691 106558.200093845
28 -8.11614527938351E-6 -86246.8400568131 -24641.2400162323 123211.200081161
29 -7.05486284239850E-6 -99221.340049384 -28348.2400141097 141746.200070549
30 -6.16142944895955E-6 -113609.00004313 -32459.0000123228 162300.000061614
31 -5.40499053405699E-6 -129508.940037835 -37001.840010810 185014.20005405
32 -4.76110670849903E-6 -147023.640033328 -42006.0400095222 210035.200047611
33 -4.21027458483498E-6 -166258.940029472 -47501.8400084205 237514.200042103
34 -3.73682023443470E-6 -187324.040026158 -53520.4400074736 267607.200037368
35 -3.32806389840558E-6 -210331.500023296 -60094.0000066561 300475.000033281
36 -2.97368408501804E-6 -235397.240020816 -67255.6400059474 336283.200029737
37 -2.66522957466352E-6 -262640.540018657 -75039.4400053305 375202.200026652
38 -2.39574209533138E-6 -292184.04001677 -83480.4400047915 417407.200023957
39 -2.15946248375913E-6 -324153.740015116 -92614.640004319 463078.200021595
40 -1.95160031217110E-6 -358679.000013661 -102479.000003903 512400.000019516

Sign in to comment.


Torsten
Torsten on 10 Feb 2023
Edited: Torsten on 10 Feb 2023
Look whether this might help.
syms T
x = sym('x',[1 4]);
eqn1 = 10.*(T-x(1))+0.2.*(T.^4-x(1).^4) - x(4) == 0;
eqn2 = -1.*(1/x(1)-1/x(2)) - x(4) == 0;
eqn3 = 2.*(x(3)-x(2))-x(4) == 0;
eqn4 = 5.*(1-x(3))-x(4) == 0;
solve([eqn1,eqn2,eqn3,eqn4],x)
Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.
ans = struct with fields:
x1: ((21*T^4)/100 + (21*T)/2 + 8)*root(z^9 - z^8*(10*T + T^4/5 + 40/7) + z^7*((400*T)/7 + (8*T^4)/7 + 390/49) - z^6*((8800*T)/49 + (176*T^4)/49 - 7200/343) + z^5*((124000*T)/343 + (2480*T^4)/343 - 273500/2401) - z^4*((1234000*T)/2401 + (24680*T^4)… x2: 1 - (7*root(z^9 - z^8*(10*T + T^4/5 + 40/7) + z^7*((400*T)/7 + (8*T^4)/7 + 390/49) - z^6*((8800*T)/49 + (176*T^4)/49 - 7200/343) + z^5*((124000*T)/343 + (2480*T^4)/343 - 273500/2401) - z^4*((1234000*T)/2401 + (24680*T^4)/2401 - 2987401/12005) … x3: 1 - root(z^9 - z^8*(10*T + T^4/5 + 40/7) + z^7*((400*T)/7 + (8*T^4)/7 + 390/49) - z^6*((8800*T)/49 + (176*T^4)/49 - 7200/343) + z^5*((124000*T)/343 + (2480*T^4)/343 - 273500/2401) - z^4*((1234000*T)/2401 + (24680*T^4)/2401 - 2987401/12005) + z… x4: root(z^9 - z^8*(10*T + T^4/5 + 40/7) + z^7*((400*T)/7 + (8*T^4)/7 + 390/49) - z^6*((8800*T)/49 + (176*T^4)/49 - 7200/343) + z^5*((124000*T)/343 + (2480*T^4)/343 - 273500/2401) - z^4*((1234000*T)/2401 + (24680*T^4)/2401 - 2987401/12005) + z^3*(…
It shows that x(1)-x(4) are related to the roots of a polynomial of degree 9 in T.

Categories

Find more on Loops and Conditional Statements 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!