数值分析 - 编程作业(5)

第八章非线性方程及非线性方程组的解法

1. 求下列方程的非零根

x?513?0.6651x?f(x)?ln???0 ?513?0.6651x1400?0.0918??

分析：本题拟采用牛顿法对该非线性方程进行求解，而本函数f(x)为奇函数，为了使得牛顿迭代公式初始值比较靠近该方程的非零根，应该确定该方程非零根的大概范围。故首先用Matlab画图指令画出该函数的曲线图。

因注意到513-0.6651x>0，则可知|x|<771.3126.故有以下绘图代码：

>> x=-770:770;

>> y=log((513+0.6651*x)./(513-0.6651*x))-x/(1400*0.0918); >> plot(x,y) >> grid

函数图如下：

3210-1-2-3-800-600-400-2000200400600800

由上图可以很明显看出，原方程有一个零根，两个关于原点对称的非零实根。由图估测原方程的正实根的x坐标值在770附近。并且经验证f(770)=1.0771>0，而f(760)= -1.0055<0。所以正实根落在区间[760,770]上。

现取初始值x0=765。误差限errorlim=10?5，最大迭代次数N=100

M文件代码如下：

%子函数非线性函数f3

function y=f3(x)

y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end

%子函数非线性函数一阶导数df3 function y=df3() syms x1

y=log((513+0.6651*x1)/(513-0.6651*x1))-x1/(1400*0.0918); y=diff(y); end

%主程序----牛顿法求解非线性方程 clear clc

format long x0=765; N=100;

errorlim=10^(-5);

x=x0-f3(x0)/subs(df3(),x0); n=1; while n

x=x0-f3(x0)/subs(df3(),x0); if abs(x-x0)>errorlim n=n+1; else break; end

x0=x; end

disp('非线性方程求解完毕！') disp(['迭代次数: n=',num2str(n)])

disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)])

运行结果：

非线性方程求解完毕！ 迭代次数: n=5

所求非零根: 正根x1=767.3861 负根x2=-767.3861

结果验证： >> f3(x) ans =

-4.440892098500626e-015 >> f3(-x) ans =

4.440892098500626e-015

第九章常微分方程数值解法

1、 设常微分方程初值问题

?y???y?2cosx,0?x?? ?y(0)?1?其精确解为y?cosx?sinx。选取步长h使四阶Adams预测-校正算法和经典RK法均稳定，分别用这两种方法求解微分方程，将数值解和精确解进行比

较，输出结果。其中多步法需要的初值由经典RK法提供。

解：

M文件代码：

%子函数----四阶RK算法

function yrk4=rk4(a,b,y0,n) h=(b-a)/n; x=a:h:b;

yrk4=zeros(n+1,1); yrk4(1)=y0;

%round—控制RK-4算法计算步数 if round~=n n=round; end

for i=1:n

k1=fxy(x(i),yrk4(i));

k2=fxy(x(i)+h/2,yrk4(i)+h*k1/2); k3=fxy(x(i)+h/2,yrk4(i)+h*k2/2); k4=fxy(x(i)+h,yrk4(i)+h*k3);

yrk4(i+1)=yrk4(i)+h*(k1+2*k2+2*k3+k4)/6; end

%子函数----fxy函数 function f=fxy(x,y) f=-y+2*cos(x);

%主程序----经典RK法求解常微分方程 clear clc a=0; b=pi;

N=50; n=1; y0=1;

yrk4=rk4(a,b,y0,N,N); %绘图---精确解

y=dsolve('Dy=-y+2*cos(x)','y(0)=1','x'); x1=0:(pi/N):pi; y1=subs(y,x1); plot(x1,y1,'-b') hold on

%绘图---数值解 plot(x1,yrk4,'xr') grid on

legend('精确解','数值解') title('RK-4数值解与精确解')

运行结果：

RK-4数值解与精确解1.5精确解数值解1 0.50-0.5-1 00.511.522.533.5

%子函数----Adams预测校正算法·¨ function yadms=adms(a,b,y,n) h=(b-a)/n; x=a:h:b;

yadms=[y;zeros(n-3,1)];