數(shù)值分析實(shí)驗(yàn)

Chapter5

SolutionofNonlinearEquationsMainpointsContentsExperiment5實(shí)驗(yàn)?zāi)康模悍蔷€性方程求解實(shí)驗(yàn)內(nèi)容：1.二分法的Matlab實(shí)現(xiàn)；

2.牛頓法的Matlab實(shí)現(xiàn)；

3.牛頓下山法、割線法、艾特金加速法、重根迭代法、非線性方程組牛頓法中任選其一。

實(shí)驗(yàn)要求：1.每種算法要求達(dá)到給定的精度，輸出近似解結(jié)果及所需迭代次數(shù)；

2.P.239、171，或自選題目；

3.每個(gè)算法至少實(shí)驗(yàn)一個(gè)題目。

求解方程（組）的solve命令MATLAB命令:>>[x1,x2,…,xn]=solve('F1=0','F2=0',...'Fn=0')Example1solvethesystem:>>[x,y,z]=solve('x^2+4*x*y+z=0','x+3*y*z=3','z^2+5*sin(z)-y=0')x=2.7235770787502765061989585805873y=-.66823751738808258767593643769513z=-.13788656181720556936532298380376%二分法的MATLAB主程序function[k,x,wuca,yx]=erfen(a,b,abtol)a(1)=a;b(1)=b;ya=fun(a(1));yb=fun(b(1));%程序中調(diào)用的fun.m為函數(shù)

ifya*yb>0,disp('注意：ya*yb>0,請(qǐng)重新調(diào)整區(qū)間端點(diǎn)a和b.'),returnendmax1=-1+ceil((log(b-a)-log(abtol))/log(2));%ceil(u)是大于u的最小取整數(shù)fork=1:max1+1a;ya=fun(a);b;yb=fun(b);x=(a+b)/2;yx=fun(x);wuca=abs(b-a)/2;k=k-1;[k,a,b,x,wuca,ya,yb,yx]ifyx==0a=x;b=x;elseifyb*yx>0b=x;yb=yx;elsea=x;ya=yx;endifb-a<abtol,return,endendk=max1;x;wuca;yx=fun(x);例2確定方程x3-x+4=0的實(shí)根的分布情況，并用二分法求在開區(qū)間(-2,-1)內(nèi)的實(shí)根的近似值，要求精度為0.001.解（1）先用兩種方法確定方程x3-x+4=0的實(shí)根的分布情況。方法1作圖法.在MATLAB工作窗口輸入如下程序>>x=-4:0.1:4;y=x.^3-x+4;plot(x,y)grid畫出函數(shù)f(x)=x3-x+4的圖像.從圖像可以看出，此曲線在(-2,-1)內(nèi)與x軸只有一個(gè)交點(diǎn)，則該方程有唯一個(gè)實(shí)根，且在(-2,-1)內(nèi).

見下頁圖。方法2試算法.在MATLAB工作窗口輸入程序>>x=-4:1:4,y=x.^3-x+4運(yùn)行后輸出結(jié)果x=-4-3-2-101234y=-56-20-2444102864由于連續(xù)函數(shù)f(x)滿足，所以此方程在(-2,-1)內(nèi)有一個(gè)實(shí)根.（2）用二分法的主程序計(jì)算.在MATLAB工作窗口輸入程序>>[k,x,wuca,yx]=erfen(-2,-1,0.001)運(yùn)行后屏幕顯示用二分法計(jì)算過程被列入表2-3，其余結(jié)果為k=9，x=-1.7959，wuca=9.7656e-004，yx=0.0037次數(shù)k左端點(diǎn)ak右端點(diǎn)bk中點(diǎn)xk函數(shù)值f(ak)函數(shù)值f(bk)函數(shù)值f(xk)0-2.0000-1.0000-1.50000.5000-2.00004.00002.12501-2.0000-1.5000-1.75000.2500-2.00002.12500.39062-2.0000-1.7500-1.87500.1250-2.00000.3906-0.71683-1.8750-1.7500-1.81250.0625-0.71680.3906-0.14184-1.8125-1.7500-1.78130.0313-0.14180.39060.12965-1.8125-1.7813-1.79690.0156-0.14180.1296-0.00486-1.7969-1.7813-1.78910.0078-0.00480.12960.06277-1.7969-1.7891-1.79300.0039-0.00480.06270.02908-1.7969-1.7930-1.79490.0020-0.00480.02900.01219-1.7969-1.7949-1.79590.0010-0.00480.01210.0037%迭代法的MATLAB主程序1function[k,piancha,xdpiancha,xk]=diedai1(x0,k)%輸入的量--x0是初始值,k是迭代次數(shù)x(1)=x0;fori=1:kx(i+1)=fun1(x(i));%程序中調(diào)用的fun1.m為函數(shù)y=φ(x)piancha=abs(x(i+1)-x(i));xdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);[(i-1)pianchaxdpianchaxk]endif(piancha>1)&(xdpiancha>0.5)&(k>3)disp('注意：此迭代序列發(fā)散，請(qǐng)重新輸入新的迭代公式')return;endif(piancha<0.001)&(xdpiancha<0.0000005)&(k>3)disp('祝賀您！此迭代序列收斂，且收斂速度較快')return;endp=[(i-1)pianchaxdpianchaxk]';例3

求方程的一個(gè)正根.在MATLAB工作窗口輸入程序>>[k,piancha,xdpiancha,xk]=diedai1(2,5)運(yùn)行后輸出用迭代公式的結(jié)果[k,piancha,xdpiancha,xk]=1.000000000000001.000000000000000.333333333333333.000000000000002.000000000000002.500000000000005.000000000000000.500000000000003.000000000000004.375000000000000.897435897435904.875000000000004.0000000000000011.757812500000001.70828603859251-6.882812500000005.0000000000000011.803741455078130.63167031671297-18.68655395507813注意：此迭代序列發(fā)散，請(qǐng)重新輸入新的迭代公式k=5，piancha=11.80374145507813，xdpiancha=0.63167031671297，xk=-18.68655395507813由以上運(yùn)行后輸出的迭代序列與根=2.31662479035540相差越來越大，即迭代序列發(fā)散，此迭代法就失敗.這時(shí)偏差piancha逐漸增大且相對(duì)誤差xdpiancha的值大于0.5.解(1)建立文件diedail.m(2)建立文件fun1.mfunctiony1=fun1(x)y1=(10-x^2)/2;用迭代公式運(yùn)行后輸出的結(jié)果[k,piancha,xdpiancha,xk]=1.000000000000000.500000000000000.200000000000002.500000000000002.000000000000000.277777777777780.125000000000002.222222222222223.000000000000000.146198830409360.061728395061732.368421052631584.000000000000000.079264426125550.034626038781162.289156626506025.000000000000000.042304047651280.018144868631152.33146067415730k=5,piancha=0.04230404765128,xdpiancha=0.01814486863115,xk=2.33146067415730可見，偏差piancha和偏差的相對(duì)誤差xdpiancha的值逐漸變小，且第5次的迭代值xk

=2.33146067415730與根=2.31662479035540接近，則迭代序列收斂，但收斂速度較慢，此迭代法較為成功.%牛頓切線法的MATLAB主程序%輸入：初始值x0;近似值xk的精度tol；f(xk)的精度tol。function[k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax)x(1)=x0;fori=1:gxmaxx(i+1)=x(i)-fnq(x(i))/(dfnq(x(i))+eps);piancha=abs(x(i+1)-x(i));xdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fnq(x(i));[(i-1)xkykpianchaxdpiancha]if(abs(yk)<ftol)&((piancha<tol)|(xdpiancha<tol))k=i-1;xk=x(i);[(i-1)xkykpianchaxdpiancha]return;endendifi>gxmaxdisp('請(qǐng)注意：迭代次數(shù)超過給定的最大值gxmax。')k=i-1;xk=x(i);[(i-1)xkykpianchaxdpiancha]return;end[(i-1),xk,yk,piancha,xdpiancha]';例3

用牛頓切線法求