我审查完公式,又从新编辑,希望你能再提点一下
我的想法是 解四个方程,a,b是参数,t(1)_t(4)为变量,t(1)_t(4)分别随a,b变化,生成三维图,
t0=[0;7;0];
a=0:5;
b=0:5;
[A,B]=meshgrid(a,b)
p1=zeros(6,6);
p2=zeros(6,6);
p3=zeros(6,6);
p4=zeros(6,6);
for i=1:length(a)
for j=1:length(b)
p(:,i)=fsolve(@(t)SH_2(t,a(i),b(j)),t0);
p1(i,j)=real(sol(1));
p2(i,j)=real(sol(2));
p3(i,j)=real(sol(3));
p4(i,j)=real(sol(4));
end
end
figure
subplot(2,2,1)
surf(A,B,p1); shading interp; colorbar
subplot(2,2,2)
surf(A,B,p2); shading interp; colorbar
subplot(2,2,3)
surf(A,B,p3); shading interp; colorbar
subplot(2,2,4)
surf(A,B,p4); shading interp; colorbar
function M = SH_2( t,a,b )
y(1)=exp(20.37788-4.60517*t(2))+exp(30.14084+log(t(1))-4.60517*t(2))+exp(38.40712+2*log(t(1))-4.60517*t(2))+exp(45.2458+3*log(t(1))-4.60517*t(2))+exp(50.31148+4*log(t(1))-4.60517*t(2))+exp(48.99901+5*log(t(1))-4.60517*t(2))+exp(2.32561-2.30259*t(2))+exp(-15.45035)+exp(-44.18661+2.30259*t(2))+exp(-74.30442+4.60517*t(2))+exp(15.24311-4.60517*t(2))+exp(15.24311-4.60517*t(2))+exp(22.03574+3*log(t(1))-2.30259*t(2))+exp(-8.77285+2*log(t(1)))+exp(-40.04195+log(t(1))+2.30259*t(2))-t(4);
y(2)=exp(log(t(1)))+exp(22.56533+log(t(1))-2.30259*t(2))+exp(30.14084+log(t(1))-4.60517*t(2))+2*exp(38.40712+2*log(t(1))-4.60517*t(2))+3*exp(45.2458+3*log(t(1))-4.60517*t(2))+4*exp(50.31148+4*log(t(1))-4.60517*t(2))+5*exp(48.99901+5*log(t(1))-4.60517*t(2))+3*exp(22.03574+3*log(t(1))-2.30259*t(2))+2*exp(-8.77285+2*log(t(1)))+exp(-40.04195+log(t(1))+2.30259*t(2))-a-2*b;
y(3)=exp(22.56533+log(t(1))-2.30259*t(2))+exp(-2.30259*t(2))-exp(-32.65066+2.30259*t(2))+2*exp(20.37788-4.60517*t(2))-2*exp(log(b))+2*exp(30.14084+log(t(1))-4.60517*t(2))+2*exp(38.40712+2*log(t(1))-4.60517*t(2))+2*exp(45.2458+3*log(t(1))-4.60517*t(2))+2*exp(50.31148+4*log(t(1))-4.60517*t(2))+2*exp(48.99901+5*log(t(1))-4.60517*t(2))+exp(2.32561-2.30259*t(2))-exp(-44.18661+2.30259*t(2))-2*exp(-74.30442+4.60517*t(2))+2*exp(15.24311-4.60517*t(2))+exp(22.03574+3*log(t(1))-2.30259*t(2))-exp(-40.04195+log(t(1))+2.30259*t(2));
M=[y(1);y(2);y(3)];
end