3-7变量非线性方程组的逆Broyden解法
matlab程序
function n_broyden
clear all %清内存
clc %清屏
format long
i=input('请输入非线性方程组个数(3-7)i=');
switch i
case 3
syms x y z
x0=input('请输入初值（三维列向量[x;y;z])=');
eps=input('请输入允许的误差精度eps=');
f1=input('请输入f1(x,y,z)=');
f2=input('请输入f2(x,y,z)=');
f3=input('请输入f3(x,y,z)=');
F=[f1;f2;f3];
J=jacobian(F,[x,y,z]); %求jacobi矩阵
Fx0=subs(F,{x,y,z},x0);
Jx0=subs(J,{x,y,z},x0);
H=inv(Jx0);
x1=x0-H*Fx0;
k=0;
while norm(x1-x0)>eps %用2范数来做循环条件
p=x1-x0;
q=subs(F,{x,y,z},x1)-subs(F,{x,y,z},x0);
H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;
x0=x1;
Fx0=subs(F,{x,y,z},x0);
x1=x1-H*Fx0;
k=k+1;
end
x1
k
case 4
syms a b c d
x0=input('请输入初值（列向量[a;b;c;d])=');
eps=input('请输入允许的误差精度eps=');
f1=input('请输入f1(a,b,c,d)=');
f2=input('请输入f2(a,b,c,d)=');
f3=input('请输入f3(a,b,c,d)=');
f4=input('请输入f4(a,b,c,d)=');
F=[f1;f2;f3;f4];
J=jacobian(F,[a,b,c,d]); %求jacobi矩阵
Fx0=subs(F,{a,b,c,d},x0);
Jx0=subs(J,{a,b,c,d},x0);
H=inv(Jx0);
x1=x0-H*Fx0;
k=0;
while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0;
q=subs(F,{a,b,c,d},x1)-subs(F,{a,b,c,d},x0);
H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;
x0=x1;
Fx0=subs(F,{a,b,c,d},x0);
x1=x1-H*Fx0;
k=k+1;
end
x1
k
case 5
syms a b c d e
x0=input('请输入初值（列向量[a;b;c;d;e])='); eps=input('请输入允许的误差精度eps=');
f1=input('请输入f1(a,b,c,d,e)=');
f2=input('请输入f2(a,b,c,d,e)=');
f3=input('请输入f3(a,b,c,d,e)=');
f4=input('请输入f4(a,b,c,d,e)=');
f5=input('请输入f5(a,b,c,d,e)=');
F=[f1;f2;f3;f4;f5];
J=jacobian(F,[a,b,c,d,e]); %求jacobi矩阵
Fx0=subs(F,{a,b,c,d,e},x0);
Jx0=subs(J,{a,b,c,d,e},x0);
H=inv(Jx0);
x1=x0-H*Fx0;
while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0;
q=subs(F,{a,b,c,d,e},x1)-subs(F,{a,b,c,d,e},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;
x0=x1;
Fx0=subs(F,{a,b,c,d,e},x0);
x1=x1-H*Fx0;
k=k+1;
end
x1
k
case 6
syms a b c d e f
x0=input('请输入初值（列向量[a;b;c;d;e;f])='); eps=input('请输入允许的误差精度eps=');
f1=input('请输入f1(a,b,c,d,e,f)=');
f2=input('请输入f2(a,b,c,d,e,f)=');
f3=input('请输入f3(a,b,c,d,e,f)=');
f4=input('请输入f4(a,b,c,d,e,f)=');
f5=input('请输入f5(a,b,c,d,e,f)=');
f6=input('请输入f6(a,b,c,d,e,f)=');
F=[f1;f2;f3;f4;f5;f6];
J=jacobian(F,[a,b,c,d,e,f]); %求jacobi矩阵
Fx0=subs(F,{a,b,c,d,e,f},x0);
Jx0=subs(J,{a,b,c,d,e,f},x0);
H=inv(Jx0);
x1=x0-H*Fx0;
k=0;
while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0;
q=subs(F,{a,b,c,d,e,f},x1)-subs(F,{a,b,c,d,e,f},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;
x0=x1;
Fx0=subs(F,{a,b,c,d,e,f},x0);
x1=x1-H*Fx0;
k=k+1;
end
k
case 7
syms a b c d e f g
format long
x0=input('请输入初值（列向量[a;b;c;d;e;f;g])=');
eps=input('请输入允许的误差精度eps=');
f1=input('请输入f1(a,b,c,d,e,f,g)=');
f2=input('请输入f2(a,b,c,d,e,f,g)=');
f3=input('请输入f3(a,b,c,d,e,f,g)=');
f4=input('请输入f4(a,b,c,d,e,f,g)=');
f5=input('请输入f5(a,b,c,d,e,f,g)=');
f6=input('请输入f6(a,b,c,d,e,f,g)=');
f7=input('请输入f7(a,b,c,d,e,f,g)=');
F=[f1;f2;f3;f4;f5;f6;f7];
J=jacobian(F,[a,b,c,d,e,f,g]); %求jacobi矩阵
Fx0=subs(F,{a,b,c,d,e,f,g},x0);
Jx0=subs(J,{a,b,c,d,e,f,g},x0);
H=inv(Jx0);
x1=x0-H*Fx0;
k=0;
while norm(x1-x0)>eps %用2范数来做循环条件
p=x1-x0;
q=subs(F,{a,b,c,d,e,f,g},x1)-subs(F,{a,b,c,d,e,f,g},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1;
x0=x1;
Fx0=subs(F,{a,b,c,d,e,f,g},x0);
x1=x1-H*Fx0;
k=k+1;
end
x1
k
end
end。