回 帖 发 新 帖 刷新版面

主题:[讨论]牛顿迭代法求非线性方程,程序老报错

程序如下:
function [z,k] = newton(y0,n,eps1)
k=1;
z=y0;
t=y0-fun0(y0)./dfun0(y0);
while abs(t-y0)>=eps1
    z=[z,t];
    y0=t;
    k=k+1;
    t=y0-fun0(y0)./dfun0(y0);
    if (k-1)>n
        error('n is full');
    end
end
end
调用了fun0和dfun0两个函数,其中
function [fun] = fun0(y)
syms y
fun=49550.*log(5000000) - 49550.*log(y - (y.*exp(y./17520000))./25 + 5000000) - (1238750.*((5000000.*log(y + 5000000))./exp(y./17520000) - (5000000.*log(y - (y.*exp(y./17520000))./25 + 5000000))./exp(y./17520000)))./y + ((991.*y)./2 + 100).*((2500.*(5000000./(exp(y./17520000).*(y + 5000000)) - (125.*log(y + 5000000))./(438.*exp(y./17520000)) + (125.*log(y - (y.*exp(y./17520000))./25 + 5000000))./(438.*exp(y./17520000)) + (5000000.*(exp(y./17520000)./25 + (y.*exp(y./17520000))./438000000 - 1))./(exp(y./17520000).*(y - (y.*exp(y./17520000))./25 + 5000000))))./y - (2500.*((5000000.*log(y + 5000000))./exp(y./17520000) - (5000000.*log(y - (y.*exp(y./17520000))./25 + 5000000))./exp(y./17520000)))./y.^2 + 2500./(exp(y./17520000).*(y + 5000000)) - (100.*(exp(y./17520000)./25 + (y.*exp(y./17520000))./438000000 - 1))./(y - (y.*exp(y./17520000))./25 + 5000000) - log(y + 5000000)./(7008.*exp(y./17520000)) + log(y - (y.*exp(y./17520000))./25 + 5000000)./(7008.*exp(y./17520000)) + (2500.*(exp(y./17520000)./25 + (y.*exp(y./17520000))./438000000 - 1))./(exp(y./17520000).*(y - (y.*exp(y./17520000))./25 + 5000000))) - (1238750.*log(y + 5000000))./exp(y./17520000) + (1238750.*log(y - (y.*exp(y./17520000))./25 + 5000000))./exp(y./17520000) + 49550;
end
但是程序运行后出现:
 [z,k] = newton(2000,40,0.0001)
??? Undefined function or method 'ge' for input arguments of type 'sym'.

Error in ==> newton at 5
while abs(t-y0)>=eps1
不知这到底是哪出现问题了,请大家帮下忙,急用

回复列表 (共4个回复)

沙发

通过验算,作了修改,请再次查错。
function [z,k] = newton(y0,n,eps1)
k=1;
z=y0;
t=y0-fun0(y0)./dfun0(y0);
while abs(t-y0)>=eps1
    z=[z,t];
    y0=t;
    k=k+1;
    t=y0-fun0(y0)./dfun0(y0);%%%%%[color=FF0000]注意此处,“dfun0(y0)”计算机不会求导,因此需要将求导表%达式列出,如下程序中通过符号求导给出。[/color]
    if (k-1)>n
        error('n is full');
    end
end
end
%调用了fun0和dfun0两个函数,其中
function [fun] = fun0(y)
%syms y      %%%注意此处,符号型的不可能求值,正如你的提示错误
fun=49550.*log(5000000) - 49550.*log(y - (y.*exp(y./17520000))./25 + 5000000) - (1238750.*((5000000.*log(y + 5000000))./exp(y./17520000) - (5000000.*log(y - (y.*exp(y./17520000))./25 + 5000000))./exp(y./17520000)))./y + ((991.*y)./2 + 100).*((2500.*(5000000./(exp(y./17520000).*(y + 5000000)) - (125.*log(y + 5000000))./(438.*exp(y./17520000)) + (125.*log(y - (y.*exp(y./17520000))./25 + 5000000))./(438.*exp(y./17520000)) + (5000000.*(exp(y./17520000)./25 + (y.*exp(y./17520000))./438000000 - 1))./(exp(y./17520000).*(y - (y.*exp(y./17520000))./25 + 5000000))))./y - (2500.*((5000000.*log(y + 5000000))./exp(y./17520000) - (5000000.*log(y - (y.*exp(y./17520000))./25 + 5000000))./exp(y./17520000)))./y.^2 + 2500./(exp(y./17520000).*(y + 5000000)) - (100.*(exp(y./17520000)./25 + (y.*exp(y./17520000))./438000000 - 1))./(y - (y.*exp(y./17520000))./25 + 5000000) - log(y + 5000000)./(7008.*exp(y./17520000)) + log(y - (y.*exp(y./17520000))./25 + 5000000)./(7008.*exp(y./17520000)) + (2500.*(exp(y./17520000)./25 + (y.*exp(y./17520000))./438000000 - 1))./(exp(y./17520000).*(y - (y.*exp(y./17520000))./25 + 5000000))) - (1238750.*log(y + 5000000))./exp(y./17520000) + (1238750.*log(y - (y.*exp(y./17520000))./25 + 5000000))./exp(y./17520000) + 49550;
end
function [fun] = dfun0(y)

fun=-49550*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)-(6193750000000/exp(1/17520000*y)/(y+5000000)-77421875/219*log(y+5000000)/exp(1/17520000*y)-6193750000000*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)+77421875/219*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y))/y+(6193750000000*log(y+5000000)/exp(1/17520000*y)-6193750000000*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y))/y^2+991/2*(12500000000/exp(1/17520000*y)/(y+5000000)-156250/219*log(y+5000000)/exp(1/17520000*y)+156250/219*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)+2500*(200000*exp(1/17520000*y)+5/438*y*exp(1/17520000*y)-5000000)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000))/y-991/2*(12500000000*log(y+5000000)/exp(1/17520000*y)-12500000000*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y))/y^2-991/2*(4*exp(1/17520000*y)+1/4380000*y*exp(1/17520000*y)-100)/(y-1/25*y*exp(1/17520000*y)+5000000)+991/2*(100*exp(1/17520000*y)+1/175200*y*exp(1/17520000*y)-2500)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)+(991/2*y+100)*((-312500/219/exp(1/17520000*y)/(y+5000000)-12500000000/exp(1/17520000*y)/(y+5000000)^2+125/3069504*log(y+5000000)/exp(1/17520000*y)+156250/219*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)-125/3069504*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)+2500*(5/219*exp(1/17520000*y)+1/1534752000*y*exp(1/17520000*y))/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)-1/7008*(200000*exp(1/17520000*y)+5/438*y*exp(1/17520000*y)-5000000)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)-2500*(200000*exp(1/17520000*y)+5/438*y*exp(1/17520000*y)-5000000)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)^2*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y)))/y-(12500000000/exp(1/17520000*y)/(y+5000000)-156250/219*log(y+5000000)/exp(1/17520000*y)+156250/219*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)+2500*(200000*exp(1/17520000*y)+5/438*y*exp(1/17520000*y)-5000000)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000))/y^2-(12500000000/exp(1/17520000*y)/(y+5000000)-156250/219*log(y+5000000)/exp(1/17520000*y)-12500000000*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)+156250/219*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y))/y^2+2*(12500000000*log(y+5000000)/exp(1/17520000*y)-12500000000*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y))/y^3-1/3504/exp(1/17520000*y)/(y+5000000)-2500/exp(1/17520000*y)/(y+5000000)^2-(1/2190000*exp(1/17520000*y)+1/76737600000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)+(4*exp(1/17520000*y)+1/4380000*y*exp(1/17520000*y)-100)/(y-1/25*y*exp(1/17520000*y)+5000000)^2*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))+1/122780160000*log(y+5000000)/exp(1/17520000*y)+1/7008*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)-1/122780160000*log(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y)+(1/87600*exp(1/17520000*y)+1/3069504000000*y*exp(1/17520000*y))/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)-1/17520000*(100*exp(1/17520000*y)+1/175200*y*exp(1/17520000*y)-2500)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)-(100*exp(1/17520000*y)+1/175200*y*exp(1/17520000*y)-2500)/exp(1/17520000*y)/(y-1/25*y*exp(1/17520000*y)+5000000)^2*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y)))+1238750*(1-1/25*exp(1/17520000*y)-1/438000000*y*exp(1/17520000*y))/(y-1/25*y*exp(1/17520000*y)+5000000)/exp(1/17520000*y);
end

板凳

改后运行结果: 


n is full

3 楼

当取参数是有结果:
[z,k] = newton(2000,4000,0.01)  %%参考验算值

请仔细修改n 和 精度。

4 楼


我照你说的错误修改了,呵呵,果真可以算出来了,非常感谢!!!

我来回复

您尚未登录,请登录后再回复。点此登录或注册