最優(yōu)化算法課程設(shè)計(jì)--解方程組及非線性方程組

1、<p><b>　　課程設(shè)計(jì)任務(wù)書(shū)</b></p><p>　　2012—2013 學(xué)年第1學(xué)期</p><p>　　課程名稱： 最優(yōu)化算法 </p><p>　　設(shè)計(jì)題目：

2、解方程組及非線性方程組 </p><p>　　完成期限：自 2012 年 12 月 17 日至 2012 年 12 月 21 日共 1 周</p><p>　　指導(dǎo)教師（簽字）： 年 月 日</p><

3、p>　　系（教研室）主任（簽字）： 年 月 日</p><p><b>　　一、問(wèn)題重述</b></p><p>　?。?） 應(yīng)用GMRES方法，用求解方程組</p><p>　?。?）分別用不動(dòng)點(diǎn)迭代方法和Newton法求解，初始迭代點(diǎn)為（0，0）</p><

4、p><b>　　二、問(wèn)題分析</b></p><p><b>　　問(wèn)題1</b></p><p>　　需用GMRES方法求解線性方程組，首先建立函數(shù)gmres函數(shù)，儲(chǔ)存如下將gmres.m保存，再新建文件并調(diào)用gmres函數(shù)，賦初始值，則可輸出函數(shù)值。</p><p><b>　　問(wèn)題2</b>

5、;</p><p>　　需用不動(dòng)點(diǎn)迭代方法和Newton法求解非線性方程組，首先建立函數(shù)fun儲(chǔ)存，如下將fun.m保存；其次建立函數(shù)dfun用來(lái)求方程組的雅克比矩陣將dfun.m保存；再編程牛頓法求解非線性方程組將newton.m保存；然后將迭代過(guò)程寫入txt文檔文件名為iteration.txt；最終在調(diào)用matlab中newton函數(shù)，并賦初始值，則可輸出函數(shù)值。</p><p>&

6、lt;b>　　三、求解方法</b></p><p><b>　　問(wèn)題1</b></p><p><b>　　clear</b></p><p><b>　　clc</b></p><p><b>　　global A;</b></p

7、><p>　　format long</p><p>　　A=[10 -1 0; -1 10 -2; 0 -2 10];</p><p>　　b=[9,7,8]'; </p><p>　　x0=[0,0,0]';</p><p>　　params=[1e-14,30]; </p><p

8、>　　fprintf (' Error_rk \n')</p><p>　　[x, error, total_iters]=gmresb(x0,b,'atv', params)</p><p>　　%fprintf('x的近似解:\n')</p><p>　　%fprintf('

9、 %10.5e\n',x)</p><p>　　%fprintf ('迭代次數(shù):\n')</p><p>　　%fprintf (' % d \n', total_iters)</p><p>　　%fprintf ('解的誤差:\n')</p><p

10、>　　%fprintf (' %11 .5e \n', error)</p><p>　　fprintf (' \n')</p><p><b>　　問(wèn)題二</b></p><p>　　function df=dfun(x);</p><p>　　%用來(lái)求解方程組的雅克比矩

11、陣儲(chǔ)存在dfun中</p><p><b>　　f=fun(x);</b></p><p>　　df=[diff(f,'x1');diff(f,'x2')];</p><p>　　df=conj(df');</p><p>　　newton([0 0],0.00001,20)<

12、;/p><p><b>　　四、求解過(guò)程</b></p><p><b>　　問(wèn)題1</b></p><p>　　function [x,error, total_iters] = gmresb(x0, b, atv, params)</p><p>　　n=length(b);</p>

13、<p>　　errtol=params(1);</p><p>　　kmax=params(2);</p><p><b>　　x=x0;</b></p><p>　　h=zeros(kmax);</p><p>　　v=zeros(n,kmax);</p><p>　　if norm(

14、x) ~=0</p><p>　　r = b-feval(atv,x);</p><p><b>　　else</b></p><p><b>　　r = b;</b></p><p><b>　　end</b></p><p>　　rho=norm(r

15、);</p><p>　　errtol=errtol*norm(b);</p><p><b>　　error=[];</b></p><p>　　error=[error,rho];</p><p>　　total_iters=0;</p><p>　　if(rho < errtol) &

16、lt;/p><p><b>　　return</b></p><p><b>　　end</b></p><p>　　v(:,1)=r/rho;</p><p><b>　　beta=rho;</b></p><p><b>　　k=0;</

17、b></p><p>　　while((rho > errtol) & (k < kmax))</p><p><b>　　k=k+1;</b></p><p>　　v(:,k+1)=feval(atv,v(:,k));</p><p><b>　　%</b></p&

18、gt;<p>　　% Modified Gram-Schmidt</p><p><b>　　%</b></p><p><b>　　for j=1:k</b></p><p>　　h(j,k)=v(:,k+1)'*v(:,j);</p><p>　　v(:,k+1)=v(:,

19、k+1)-h(j,k)*v(:,j);</p><p><b>　　end</b></p><p>　　h(k+1,k)=norm(v(:,k+1));</p><p><b>　　%</b></p><p>　　% Watch out for happy breakdown.</p>

20、;<p><b>　　%</b></p><p>　　if(h(k+1,k) ~= 0)</p><p>　　v(:,k+1)=v(:,k+1)/h(k+1,k);</p><p><b>　　end</b></p><p>　　% Formulate and solve the le

21、ast squares problem.</p><p>　　% Update the residual approximation.</p><p><b>　　%</b></p><p>　　y=h(1:k+1,1:k)\(beta*eye(k+1,1));</p><p>　　rho=norm(beta*eye(k

22、+1,1) - h(1:k+1,1:k)*y);</p><p>　　error=[error,rho];</p><p><b>　　end</b></p><p>　　% At this point either k > kmax or rho < errtol.</p><p>　　% It's

23、 time to compute x and leave.</p><p><b>　　%</b></p><p>　　total_iters=k;</p><p>　　x = x + v(1:n,1:k)*y;</p><p>　　function [ax] = atv(x)</p><p>&

24、lt;b>　　global A</b></p><p><b>　　ax=A*x;</b></p><p><b>　　問(wèn)題二</b></p><p>　　function x=newton(x0,eps,N);</p><p><b>　　con=0;</b>

25、;</p><p>　　%其中x0為迭代初值eps為精度要求N為最大迭代步數(shù)con用來(lái)記錄結(jié)果是否收斂</p><p>　　for i=1:N;</p><p>　　f=subs(fun(x0),{'x1' 'x2'},{x0(1) x0(2)});</p><p>　　df=subs(dfun(x0),{&#

26、39;x1' 'x2'},{x0(1) x0(2)});</p><p>　　x=x0-f/df;</p><p>　　for j=1:length(x0);</p><p>　　il(i,j)=x(j);</p><p><b>　　end</b></p><p>　　i

27、f norm(x-x0)<eps</p><p><b>　　con=1;</b></p><p><b>　　break;</b></p><p><b>　　end</b></p><p><b>　　x0=x;</b></p>&

28、lt;p><b>　　end</b></p><p><b>　　五、結(jié)果</b></p><p><b>　　問(wèn)題一</b></p><p>　　Error_rk </p><p>　　x =1.00000000000000</p><p>

29、　　1.00000000000000</p><p>　　1.00000000000000</p><p>　　error=13.92838827718412 1.8675025938955 0.19666640169080 0.00000000000000</p><p>　　total_iters = 3</p><p>&l

30、t;b>　　>></b></p><p><b>　　問(wèn)題二</b></p><p><b>　　ans =</b></p><p>　　1.54634284935363 1.39117624667800</p><p><b>　　六、附件</b&

31、gt;</p><p><b>　　問(wèn)題一</b></p><p>　　function [x,error, total_iters] = gmresb(x0, b, atv, params)</p><p>　　n=length(b);</p><p>　　errtol=params(1);</p>&l

32、t;p>　　kmax=params(2);</p><p><b>　　x=x0;</b></p><p><b>　　%%</b></p><p>　　h=zeros(kmax);</p><p>　　v=zeros(n,kmax);</p><p>　　%r=b-f

33、eval(atv,x);</p><p>　　if norm(x) ~=0</p><p>　　r = b-feval(atv,x);</p><p><b>　　else</b></p><p><b>　　r = b;</b></p><p><b>　　end

34、</b></p><p>　　rho=norm(r);</p><p>　　errtol=errtol*norm(b);</p><p>　　error=[];%</p><p>　　% test for termination on entry%</p><p>　　error=[error,rho];&

35、lt;/p><p>　　total_iters=0;</p><p>　　if(rho < errtol) </p><p><b>　　return</b></p><p><b>　　end%%</b></p><p>　　v(:,1)=r/rho;</p>

36、<p><b>　　beta=rho;</b></p><p><b>　　k=0;%</b></p><p>　　% GMRES iteration</p><p><b>　　%</b></p><p>　　while((rho > errtol) &a

37、mp; (k < kmax))</p><p><b>　　k=k+1;</b></p><p>　　v(:,k+1)=feval(atv,v(:,k));%</p><p>　　% Modified Gram-Schmidt%</p><p><b>　　for j=1:k</b></

38、p><p>　　h(j,k)=v(:,k+1)'*v(:,j);</p><p>　　v(:,k+1)=v(:,k+1)-h(j,k)*v(:,j);</p><p><b>　　end</b></p><p>　　h(k+1,k)=norm(v(:,k+1));%</p><p>　　%

39、 Watch out for happy breakdown.%</p><p>　　if(h(k+1,k) ~= 0)</p><p>　　v(:,k+1)=v(:,k+1)/h(k+1,k);</p><p><b>　　end</b></p><p><b>　　%</b></p>

40、<p>　　% Formulate and solve the least squares problem.</p><p>　　% Update the residual approximation.</p><p><b>　　%</b></p><p>　　y=h(1:k+1,1:k)\(beta*eye(k+1,1));&

41、lt;/p><p>　　rho=norm(beta*eye(k+1,1) - h(1:k+1,1:k)*y);</p><p>　　error=[error,rho];</p><p><b>　　end</b></p><p>　　%% At this point either k > kmax or rho <

42、; errtol.</p><p>　　% It's time to compute x and leave.</p><p><b>　　%</b></p><p>　　total_iters=k;</p><p>　　x = x + v(1:n,1:k)*y;</p><p>　　fu

43、nction [ax] = atv(x)</p><p><b>　　global A</b></p><p><b>　　ax=A*x;</b></p><p><b>　　clear</b></p><p><b>　　clc</b></p>

44、;<p><b>　　global A;</b></p><p>　　format long</p><p>　　A=[10 -1 0; -1 10 -2; 0 -2 10];</p><p>　　b=[9,7,8]'; </p><p>　　x0=[0,0,0]';</p>

45、<p>　　params=[1e-14,30]; </p><p>　　fprintf (' Error_rk \n')</p><p>　　[x, error, total_iters]=gmresb(x0,b,'atv', params)</p><p>　　%fprintf('x的近似解

46、:\n')</p><p>　　%fprintf(' %10.5e\n',x)</p><p>　　%fprintf ('迭代次數(shù):\n')</p><p>　　%fprintf (' % d \n', total_iters)</p><p>

47、;　　%fprintf ('解的誤差:\n')</p><p>　　%fprintf (' %11 .5e \n', error)</p><p>　　fprintf (' \n')</p><p><b>　　問(wèn)題二</b></p><p>　　function

48、 f=fun(x);</p><p>　　syms x1 x2 </p><p>　　f1=x1^2-x2-1;</p><p>　　f2=(x1-2)^2+(x2-0.5)^2-1;</p><p>　　f=[f1 f2];</p><p><b>　　建立函數(shù)dfun</b></p>

49、;<p>　　用來(lái)求方程組的雅克比矩陣將dfun.m保存到工作路徑中:</p><p>　　function df=dfun(x);</p><p>　　%用來(lái)求解方程組的雅克比矩陣儲(chǔ)存在dfun中</p><p><b>　　f=fun(x);</b></p><p>　　df=[diff(f,'

50、x1');diff(f,'x2')];</p><p>　　df=conj(df');</p><p>　　編程牛頓法求解非線性方程組將newton.m保存到工作路徑中:</p><p>　　function x=newton(x0,eps,N);</p><p><b>　　con=0;</b

51、></p><p>　　%其中x0為迭代初值eps為精度要求N為最大迭代步數(shù)con用來(lái)記錄結(jié)果是否收斂</p><p>　　for i=1:N;</p><p>　　f=subs(fun(x0),{'x1' 'x2'},{x0(1) x0(2)});</p><p>　　df=subs(dfun(x0),

52、{'x1' 'x2'},{x0(1) x0(2)});</p><p>　　x=x0-f/df;</p><p>　　for j=1:length(x0);</p><p>　　il(i,j)=x(j);</p><p><b>　　end</b></p><p>

53、　　if norm(x-x0)<eps</p><p><b>　　con=1;</b></p><p><b>　　break;</b></p><p><b>　　end</b></p><p><b>　　x0=x;</b></p>

54、<p><b>　　end</b></p><p>　　將迭代過(guò)程寫入txt文檔文件名為iteration.txt</p><p>　　fid=fopen('iteration.txt','w');</p><p>　　fprintf(fid,'iteration');</p&

55、gt;<p>　　for j=1:length(x0)</p><p>　　fprintf(fid,' x%d',j);</p><p><b>　　end</b></p><p><b>　　for j=1:i</b></p><p>　　fprint

56、f(fid,'\n%6d ',j);</p><p>　　for k=1:length(x0)</p><p>　　fprintf(fid,' %10.6f',il(j,k));</p><p><b>　　end</b></p><p><b>　　end</b&

57、gt;</p><p><b>　　if con==1</b></p><p>　　fprintf(fid,'\n計(jì)算結(jié)果收斂！');</p><p><b>　　end</b></p><p><b>　　if con==0</b></p>&l