当前位置: 首页 > 程序语言 >MATLAB求解非线性方程的NEWTON迭代程序

MATLAB求解非线性方程的NEWTON迭代程序

作者 西入秦关
来源: 小木虫 350 7 举报帖子
+关注

用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左右

  • yanze

    for循环下面第一个判断条件abs(subs(diff(fun,'x'),x,t1))
    这个表达式是符号表达 ,每一次循环都加长了许多,后来就越来越慢,循环到第六次基本就不动了。
    你这个方法根本不是数值求解。不要出现符号表达式!

  • yanze

    参考这个吧:https://jingyan.baidu.com/article/ab69b270f5a6702ca7189f38.html

    如果对求解方法没有要求的话直接用内置的fsolve函数更好,选择马夸特算法效果不错的

  • 西入秦关

    引用回帖:
    4楼: Originally posted by yanze at 2017-12-26 20:45:57
    参考这个吧: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)

  • ecnufsy

    最好不要在数值运算上使用符号运算,除非是前处理只做一次的运算。
    里面fun函数的导数有两种选择,一是推导一下显示的导数,二是采用差分格式,如下,
    [latex]y = f(x)[/latex]
    [latex]\frac{dy}{dx} = \frac{f(x+\epsilon\Delta x) - f(x-\epsilon\Delta x)}{2\epsilon}[/latex]

  • yanze

    引用回帖:
    6楼: Originally posted by ecnufsy at 2017-12-26 22:46:24
    最好不要在数值运算上使用符号运算,除非是前处理只做一次的运算。
    里面fun函数的导数有两种选择,一是推导一下显示的导数,二是采用差分格式,如下,
    y = f(x)
    \frac{dy}{dx} = \frac{f(x+\epsilon\Delta x) - f ...

    接着回答,他的这个问题并不像看起来那么简单,方程里面有个系数是级数求和
    我用了fsolve函数的马夸特算法,优点就是收敛快。
    程序在附件里。
    结果是:  k=1      35.4497

  • 西入秦关

    引用回帖:
    6楼: Originally posted by ecnufsy at 2017-12-26 22:46:24
    最好不要在数值运算上使用符号运算,除非是前处理只做一次的运算。
    里面fun函数的导数有两种选择,一是推导一下显示的导数,二是采用差分格式,如下,
    y = f(x)
    \frac{dy}{dx} = \frac{f(x+\epsilon\Delta x) - f ...

    谢谢回帖

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓