MATLAB求解非线性方程的NEWTON迭代程序
用newton迭代求一直处于正忙
function x1=newton(t1,esp,mk)
k=1;
C=0;
p0=100;
E=26000;
v=0.3;
G=E/2/(1+v);
M=E/(1-v*v*(2-k));
phi=30*pi/180;
psi=0.5*pi/180;
Y=2*C*cos(phi)/(1-sin(phi));
mm=2*G/p0;
Rps=(1+sin(phi))/(1-sin(phi));
beta=(1+sin(psi))/(1-sin(psi));
delta=(Y+(Rps-1)*p0)/(2*(k+Rps)*G);
xi=(1/M)*(beta-(k*v/(1-v*(2-k)))+1/Rps/M*(k*(1-2*v)+2*v-k*beta*v/(1-v*(2-k))));
q=((1+k)*Rps*(Y+(Rps-1)*p0))/(Rps-1)/(k+Rps);
s=-(xi*(Rps-1)*q^k)/(Rps*beta);
gamma=(Rps*beta*s)/((k*Rps)-(k*beta*(Rps-1))+Rps*beta);
lambda=(Rps-1)/Rps
n0=0;
n=100;
syms t An
An=1/factorial(t)*((xi*q/beta)^t)*(Rps*beta*s)/(k*Rps-k*beta*(Rps-1)*(1+t)+Rps*beta)
AA=double(symsum(An,t,n0,inf))
A1=xi*q/beta
syms x;
fun=(exp(-A1*(x)^(k*lambda))*(AA*(x)^(k*lambda*(n+1)-1)+(delta*(1+k)*exp(A1)-AA)*(x)^(k/beta)))-(x)^(-1);
for kk=1:mk
if abs(subs(diff(fun,'x'),x,t1))<esp
x1=t1;
break;
else
if subs(diff(fun,'x',2),x,t1)==0
break;
disp('解题失败!')
else
t0=t1;
t1=t0-subs(fun,x,t0)/subs(diff(fun,'x'),x,t0);
if abs(t1-t0)<esp
x1=t1;
break;
end
end
end
end
x1=vpa(x1,8);
%求解
x1=newton(1,1e-4,100)
原计算公式
结果应该35左右
for循环下面第一个判断条件abs(subs(diff(fun,'x'),x,t1))
这个表达式是符号表达 ,每一次循环都加长了许多,后来就越来越慢,循环到第六次基本就不动了。
你这个方法根本不是数值求解。不要出现符号表达式!
参考这个吧:https://jingyan.baidu.com/article/ab69b270f5a6702ca7189f38.html
如果对求解方法没有要求的话直接用内置的fsolve函数更好,选择马夸特算法效果不错的
谢谢您的回帖,但是,我刚用fsolve函数,还是没解决,能再指教一下吗?谢谢了
function F=y(x)
k=1;
C=0;
p0=100;
E=26000;
v=0.3;
G=E/2/(1+v);
M=E/(1-v*v*(2-k));
phi=30*pi/180;
psi=30*pi/180;
Y=2*C*cos(phi)/(1-sin(phi));
mm=2*G/p0;
Rps=(1+sin(phi))/(1-sin(phi));
beta=(1+sin(psi))/(1-sin(psi));
delta=(Y+(Rps-1)*p0)/(2*(k+Rps)*G);
xi=(1/M)*(beta-(k*v/(1-v*(2-k))))+1/Rps/M*(k*(1-2*v)+2*v-k*beta*v/(1-v*(2-k)));
q=((1+k)*Rps*(Y+(Rps-1)*p0))/(Rps-1)/(k+Rps);
s=-(xi*(Rps-1)*q*k)/(Rps*beta);
%gamma=(Rps*beta*s)/((k*Rps)-(k*beta*(Rps-1))+Rps*beta);
lambda=(Rps-1)/Rps
n0=0;
n=100;
syms t An
An=1/factorial(t)*((xi*q/beta)^t)*(Rps*beta*s)/(k*Rps-k*beta*(Rps-1)*(1+t)+Rps*beta)
AA=double(symsum(An,t,n0,inf))
A1=xi*q/beta
F=(exp(-A1*(x)^(k*lambda))*(AA*(x)^(k*lambda*(n+1)-1)+(delta*(1+k)*exp(A1)-AA)*(x)^(k/beta)))-(x)^(-1);
end
%
clc
clear all
x0=1;
option=optimset('Display','iter')
x=fsolve(@(x)y(x),x0,option)
q=y(x)
,
最好不要在数值运算上使用符号运算,除非是前处理只做一次的运算。
里面fun函数的导数有两种选择,一是推导一下显示的导数,二是采用差分格式,如下,
[latex]y = f(x)[/latex]
[latex]\frac{dy}{dx} = \frac{f(x+\epsilon\Delta x) - f(x-\epsilon\Delta x)}{2\epsilon}[/latex]
接着回答,他的这个问题并不像看起来那么简单,方程里面有个系数是级数求和
我用了fsolve函数的马夸特算法,优点就是收敛快。
程序在附件里。
结果是: k=1 35.4497
谢谢回帖