matlab利用for循环及fsolve函数解数组方程

2023-07-02 09:01| 来源: 网络整理| 查看: 265

matlab利用for循环及fsolve函数解数组方程。如Ax+B=C，A,B,C均为向量，A=[1,2,3,4,5,6…]',求出x的所有值。本文中需要求解的方程较为复杂，如图片所示。NN、MM、t、v、l均为列向量，求解T2、P、N2，其余参数均为常数。需求解的方程如图片中所示结合for循环与fsolve函数对方程进行求解。函数文件如下（表达式中方程均通过移项改写为不含除法的等式）：

function eq=myfun2(p) global t v l MM NN i S2=0.0499; m=70; g=9.8; Pa=0.1*10.^6; R=320; V20=0.0069; X2=0.75; k=1.219; T1=2949; Q=1.03; K0=0.6522; Q22=0.4; St=0.000397; P=p(1); N2=p(2); T2=p(3); eq(1)=((X2*k*T1*R.*NN(i))-(k-1)*S2.*P.*v(i))-T2.*(R.*(NN(i)+(k-1).*N2)); eq(2)=((MM(i)-N2.*t(i)).*T2.*R)-P.*(V20+S2.*l(i)); eq(3)=(Q22*K0.*P.*St)-N2.*((R.*T2).^(0.5)); end

脚本文件如下：

global t v l MM NN i t=[0 7.46798E-06 1.4936E-05 2.24039E-05 2.98719E-05 6.72118E-05 0.000104552 0.000141892 0.000179232 0.000365931 0.000552631 0.00073933]'; v=[0 -5.00678E-05 -9.97945E-05 -0.000149178 -0.000198217 -0.000438178 -0.000669247 -0.000891203 -0.001103826 -0.002019452 -0.002670596 -0.003033269]'; l=[0 -1.87164E-10 -7.46962E-10 -1.67684E-09 -2.97423E-09 -1.4883E-08 -3.55866E-08 -6.47488E-08 -1.02025E-07 -3.97527E-07 -8.39624E-07 -1.37672E-06]'; MM=[8.435E-09 5.34858E-07 1.06428E-06 1.59667E-06 2.13205E-06 4.85335E-06 7.64797E-06 1.05149E-05 1.34534E-05 2.91855E-05 4.65705E-05 6.55032E-05]'; NN=[0.0440100000000000 0.0444490931018484 0.0448872450772039 0.0453244479498637 0.0457607159982832 0.0479279889460761 0.0500719706938977 0.0521928056854323 0.0542906780878086 0.0644415463843334 0.0740439392842003 0.0831186265744276]'; TT=zeros(12,1); PP=zeros(12,1); N22=zeros(12,1); exit0=zeros(12,1); f1=zeros(12,1); f2=zeros(12,1); f3=zeros(12,1); for i=1:1:12 [p,fval,exitflag]=fsolve('myfun2',[100000,0.005,2000],optimset('Display','off')); PP(i)=p(1); N22(i)=p(2); TT(i)=p(3); exit0(i)=exitflag; f1(i)=fval(1); f2(i)=fval(2); f3(i)=fval(3); end

运行以上代码就能得到每一次循环中T2、P、N2的值，并且通过fval及exitflag的值能够判断出所求得解是否可信。初值的设定在一定程度上会影响解的结果，因此初值需要尽量设定在真解附近。也可通过fval及exitflag的值来修正初值，直至fval及exitflag为理想值为止。可能需要注意几点： 1.一定要在函数文件及脚本文件都设置t、v、l、MM、NN、i为全局变量； 2.函数文件中尽可能不写除号，可将除法方程移项改写为乘法方程； 3.for循环一定要写在脚本文件里，写在函数文件里只能得到最后一次的运行结果。本人自学小白，没有把握确定以上代码完全正确，仅为分享，轻喷。若有欠缺错误，欢迎指正并一同探讨学习。

【本文地址】

公司简介

联系我们