Question

有哪位matlab高手帮我看下下面这段程序，是利用四个方程解四个未知数

syms k k1 k2 k3
s=150,d=3800,d1=1200,d2=1500,d3=1100,b=122,b1=80,b2=60,b3=70,sl=0.96,h=6,h1=12,h2=15,h3=10,l11=0.01,l12=0.01,l13=0.02,l01=0.15,l02=0.16,l03=0.12,s1=50,s2=60,s3=40
a1=0.436183,a2=-0.120167,a3=0.9372986,a=0.33267
t=1/(1+a*k);
t1=1/(1+a*k1);
t2=1/(1+a*k2);
t3=1/(1+a*k3);
f=(1/(2*pi)^(1/2))*exp(-1/2*k^2);
g=1-1/(2*pi)^(1/2)*exp(-k^2/2)*(a1*t+a2*t^2+a3*t^3);
f1=(1/(2*pi)^(1/2))*exp(-1/2*k1^2);
g1=1-1/(2*pi)^(1/2)*exp(-k1^2/2)*(a1*t1+a2*t1^2+a3*t1^3);
f2=(1/(2*pi)^(1/2))*exp(-1/2*k2^2);
g2=1-1/(2*pi)^(1/2)*exp(-k2^2/2)*(a1*t2+a2*t2^2+a3*t2^3);
f3=(1/(2*pi)^(1/2))*exp(-1/2*k3^2);
g3=1-1/(2*pi)^(1/2)*exp(-k3^2/2)*(a1*t3+a2*t3^2+a3*t3^3);
p=k*g+f-k;
q=1+k*f-g+k^2*g-k^2*g^2-f^2-2*k*f*g;
p1=k1*g1+f1-k1;
q1=1+k1*f1-g1+k1^2*g1-k1^2*g1^2-f1^2-2*k1*f1*g1;
p2=k2*g2+f2-k2;
q2=1+k2*f2-g2+k2^2*g2-k2^2*g2^2-f2^2-2*k2*f2*g2;
p3=k3*g3+f3-k3;
q3=1+k3*f3-g3+k3^2*g3-k3^2*g3^2-f3^2-2*k3*f3*g3;
F=-s*d*(1-sl)*(g-1)/(b*(k*g+f-k)^2)+h*b*(g-1)/(2*(1-sl))+h*b
G=(1/2)*h1*k1*((l01+l11*b1)*p1+(d1)^2*(l01)^2*b1*q1)^(-1/2)*((l01+l11*b1)*(g1-1)+(d1)^2*(l01)^2*b1*(2*k1*g1-2*k1*g1^2-2*k1*f1^2-2*f1*g1-2*k1^2*f1*g1))+(1/2)*h2*k2*((l02+l12*b2)*p2+(d2)^2*(l02)^2*b2*q2)^(-1/2)*((l02+l12*b2)*(g2-1)+(d2)^2*(l02)^2*b2*(2*k2*g2-2*k2*g2^2-2*k2*f2^2-2*f2*g2-2*k2^2*f2*g2))+(1/2)*h3*k3*((l03+l13*b3)*p3+(d3)^2*(l03)^2*b3*q3)^(-1/2)*((l03+l13*b3)*(g3-1)+(d3)^2*(l03)^2*b3*(2*k3*g3-2*k3*g3^2-2*k3*f3^2-2*f3*g3-2*k3^2*f3*g3))
H=F+G
I1=-s1*d1*(1-sl)*(g1-1)/(b1*(k1*g1+f1-k1)^2)+h1*b1*(g1-1)/(2*(1-sl))+h1*b1
I2=-s2*d1*(1-sl)*(g2-1)/(b2*(k2*g2+f2-k2)^2)+h2*b2*(g2-1)/(2*(1-sl))+h2*b2
I3=-s3*d1*(1-sl)*(g3-1)/(b3*(k3*g3+f3-k3)^2)+h3*b3*(g3-1)/(2*(1-sl))+h3*b3
[k,k1,k2,k3]=solve(H,I1,I2,I3)
上面这个程序能运行，但是有点小问题，
有matlab的可以尝试运行下，帮我找下问题在哪，不胜感激

Answer 1

你把我的matlab憋死好几次。

clc;clear
syms k k1 k2 k3
s=150,d=3800,d1=1200,d2=1500,d3=1100,b=122,b1=80,b2=60,b3=70,sl=0.96,h=6,h1=12,h2=15,h3=10,l11=0.01,l12=0.01,l13=0.02,l01=0.15,l02=0.16,l03=0.12,s1=50,s2=60,s3=40
a1=0.436183,a2=-0.120167,a3=0.9372986,a=0.33267
t=1/(1+a*k);
t1=1/(1+a*k1);
t2=1/(1+a*k2);
t3=1/(1+a*k3);
f=(1/(2*pi)^(1/2))*exp(-1/2*k^2);
g=1-1/(2*pi)^(1/2)*exp(-k^2/2)*(a1*t+a2*t^2+a3*t^3);
f1=(1/(2*pi)^(1/2))*exp(-1/2*k1^2);
g1=1-1/(2*pi)^(1/2)*exp(-k1^2/2)*(a1*t1+a2*t1^2+a3*t1^3);
f2=(1/(2*pi)^(1/2))*exp(-1/2*k2^2);
g2=1-1/(2*pi)^(1/2)*exp(-k2^2/2)*(a1*t2+a2*t2^2+a3*t2^3);
f3=(1/(2*pi)^(1/2))*exp(-1/2*k3^2);
g3=1-1/(2*pi)^(1/2)*exp(-k3^2/2)*(a1*t3+a2*t3^2+a3*t3^3);
p=k*g+f-k;
q=1+k*f-g+k^2*g-k^2*g^2-f^2-2*k*f*g;
p1=k1*g1+f1-k1;
q1=1+k1*f1-g1+k1^2*g1-k1^2*g1^2-f1^2-2*k1*f1*g1;
p2=k2*g2+f2-k2;
q2=1+k2*f2-g2+k2^2*g2-k2^2*g2^2-f2^2-2*k2*f2*g2;
p3=k3*g3+f3-k3;
q3=1+k3*f3-g3+k3^2*g3-k3^2*g3^2-f3^2-2*k3*f3*g3;
F=-s*d*(1-sl)*(g-1)/(b*(k*g+f-k)^2)+h*b*(g-1)/(2*(1-sl))+h*b ;
G=(1/2)*h1*k1*((l01+l11*b1)*p1+(d1)^2*(l01)^2*b1*q1)^(-1/2)*((l01+l11*b1)*(g1-1)+(d1)^2*(l01)^2*b1*(2*k1*g1-2*k1*g1^2-2*k1*f1^2-2*f1*g1-2*k1^2*f1*g1))+(1/2)*h2*k2*((l02+l12*b2)*p2+(d2)^2*(l02)^2*b2*q2)^(-1/2)*((l02+l12*b2)*(g2-1)+(d2)^2*(l02)^2*b2*(2*k2*g2-2*k2*g2^2-2*k2*f2^2-2*f2*g2-2*k2^2*f2*g2))+(1/2)*h3*k3*((l03+l13*b3)*p3+(d3)^2*(l03)^2*b3*q3)^(-1/2)*((l03+l13*b3)*(g3-1)+(d3)^2*(l03)^2*b3*(2*k3*g3-2*k3*g3^2-2*k3*f3^2-2*f3*g3-2*k3^2*f3*g3)) ;
H=F+G ;
I1=-s1*d1*(1-sl)*(g1-1)/(b1*(k1*g1+f1-k1)^2)+h1*b1*(g1-1)/(2*(1-sl))+h1*b1 ;
I2=-s2*d1*(1-sl)*(g2-1)/(b2*(k2*g2+f2-k2)^2)+h2*b2*(g2-1)/(2*(1-sl))+h2*b2 ;
I3=-s3*d1*(1-sl)*(g3-1)/(b3*(k3*g3+f3-k3)^2)+h3*b3*(g3-1)/(2*(1-sl))+h3*b3 ;
H=vpa(H,5)
I1=vpa(I1,5)
I2=vpa(I2,5)
I3=vpa(I3,5)
k1= solve(I1)
k2= solve(I2)
k3= solve(I3)
H=subs(H);
k=solve(H)
%[k,k1,k2,k3]=solve(H,I1,I2,I3)

结果：
k1 =

1.0573072862212452544182530973701

k2 =

.95885784931743624416550059139032

k3 =

1.0172388898500092448466743483510

k =

2.0738748605418638038863473122653