数值分析实验九.docVIP

实验九 一.实验目的：数值积分、柯西问题和非线性方程的算法设计 及实现。 二.实验要求： 分别用矩形法、梯形法和西普森公式三种算法实现； 分别euler格式和预估矫正格式对ODE进行数值求解； 分别采用一般迭代格式和牛顿迭代格式对非线性方程进行数值求解； 对上述三部分的实例请选用各章节实验题一 实验正文 一．矩形法、梯形法和西普森公式程序: function y=f(x) y=log(1+x).*log(1-x); 1. 西普森公式法 function y=simpson(f,a,b,n) h=(b-a)/(2*n); y1=0; y2=0; for k=1:n x0=a+h*(2*k-1); y1=y1+limit(f,x0); end for k=1:(n-1) x0=a+h*k*2; y2=y2+limit(f,x0); end y=h*(limit(f,a)+limit(f,b)+4*y1+2*y2)/3; y=double(y); 运行结果： y=simpson(f,-1,1,100) y = -1.1016 2. 梯形法 function y=tx(f,a,b,n) h=(b-a)/n; y=0; for k=1:(n-1) x0=a+h*k; y=y+limit(f,x0); end y=h*(limit(f,a)+limit(f,b))/2+h*y; y=double(y); 运行结果： y=tx(f,-1,1,100) y = -1.1016 二．ODE的数值求解 Euler算法1： function y=euler1(a,b,f,y0,h) n=(b-a)/h; y=ones(1,n+1); y(1)=y0; x=a:h:b; for i=1:n y(i+1)=y(i)+h*feval(f,x(i),y(i)); end Euler算法2： function y=euler2(a,b,f,y0,h) n=(b-a)/h; y=ones(1,n+1); y(1)=y0; x=a:h:b; for i=1:n z=y(i)+h*feval(f,x(i),y(i)); y(i+1)=y(i)+h*( feval(f,x(i),y(i)+ feval(f,x(i+1),z))/2; end 运行结果： 1）步长h=0.2； f=inline(y-20,x,y); y1=euler1(0,10,f,10,0.2); y2=euler2(0,10,f,10,0.2); x=0:0.2:10; y=20+exp(x); plot(x,y,k,x,y1,k:,x,y2,*); 2）步长h=0.5； f=inline(y-20,x,y); y1=euler1(0,10,f,10,0.5); y2=euler2(0,10,f,10,0.5); x=0:0.5:10; y=20+exp(x); plot(x,y,k,x,y1,k:,x,y2,*); 实验结果分析：上述实验所取步长h分别为0.2和0.5，从实验结果可以看出当步长取0.5时所得数值解逼近原微分方程的精确程度更好。 三. 非线性方程进行数值求解 1.一般迭代 编制函数文件 1)function y=g(x) y=log(x.^3); 2)function y=new(x) x1=g(x); n=1; while (abs(x1-x)1.0e-6)(n=1000) x=x1;n=n+1; end x1 n 运行结果： new(0) Warning: Log of zero. In g at 2 In new at 2 x1 = -Inf n = 2 new(1) x1 = 0 n = 2 new(2) x1 = 2.0794 n = 2 2.牛顿迭代 function y=newton_1(a,n,x0,nn,eps1) x(1)=x0； b=1； i=1； while(abs(b)eps1*x(i)) i=i+1； x(i)=x(i-1)-n_f(a,n,x(i-1))/n_df(a,n,x(i-1))； b=x(i)-x(i-1)； if(inn)error(ˊnn is fullˊ)； return； end end y=x(i)； i end 程序中调用的n_f.m和n_df.m文件如下： function y=n_df(a,n,x) y=0.0; for i=1:n y=y+a(i)*(n+1-i)*x^(n-i); end function y=n_df(a,n,x) y