数值分析课程设计.doc

数值分析课程设计一、Gauss消元法1.1算法原理对n阶线性方程组Ax=b，记,，将方程组写成增广矩阵的形式=设，依次用乘矩阵的第1行加到第i行，经过k-1步消元过程，当消元过程进行到第k步时，其相应的增广矩阵为= 选择第k列的元素中绝对值最大的元素作为主元，记令 若，此时矩阵不可逆，方程的解不确定（无解或有无穷多解），则停止计算；否则，当rk时，则在其增广矩阵中交换第k行和第r行，即 使成为主元（此时变成了）。然后按Gauss消去法进行消元运算。1.2实验目的与要求1 熟悉Gauss消元法的基本原理.2 能编程实现Gauss消元法，并求解线性方程组。3运用Matlab编程，输入每个线形方程组的系数和常数，求出线形方程组的解1.3算法程序function x,det,index=Gauss(A,b)% 求线性方程组的列主元Gauss消去法，其中，% A为方程组的系数矩阵；% b为方程组的右端项；% x为方程组的解；% det为系数矩阵A的行列式的值；% index为指标变量，index=0表示计算失败；% index=1表示计算成功。n,m=size(A);nb=length(b);% 当方程组行与列的维数不相等时的，停止计算，并输出出错信息。if n=merror(The rows and columns of matrix A must be equal!);return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息。if n=nberror(The columns of A must be equal the length of b!);return;end%开始计算，先赋初值index=1;det=1;x=zeros(n,1);%给行列式的值赋1，生成n个零元素的列阵for k = 1:n-1 %k的值取1，2n-1%选主元a_max= 0;for i = k:nif abs(A(i,k)a_max %如果矩阵A中第i行第k列元素的绝对值大于a_max=abs(A(i,k);r=i; % a_max，则把此值赋给a_maxendendif a_maxkfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z; %两行中对应每个元素互换endz=b(k);b(k)=b(r);b(r)=z;det=-det; %两行互换，行列式变负号end%消元过程for i=k+1:nm=A(i,k)/A(k,k); %某列中列主元素以下中的变量系数分别除于此主元素for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k); %行列式的值乘列主元enddet=det*A(n,n); %回代过程if abs(A(n,n)A=1 1 1;1 3 -2;5 -1 4;b=6 1 19;x,det,indet=Gauss(A,b)x = 2.0000 1.4000 2.6000det =-20.0000indet =12求解线性方程组用Matlab求解A=-3 2 6;10 -7 0;5 -1 5;b=4 7 6;x,det,indet=Gauss(A,b)x =0.0000 -1.0000 1.0000det = 155indet =1二 Lagrange插值2.1算法原理构造n次多项式，其中基函数 ，显然满足 ，此时，误差 其中且依赖于，很显然，当n=1、插值节点只有两个，时2.2实验目的与要求1 熟悉Lagrange插值的基本原理；2能编程实现Lagrange插值，并求解函数多项式的值；3运用Matlab编程，根据实例中给定的函数值表求出插值多项式和函数在某一点的近似值 2.3算法程序function yi=Lagrange(x,y,xi)% Lagrange插值多项式，其中，% x为向量，全部的插值节点；% y为向量，插值节点处的函数值；% xi为标量，被估计函数的自变量：% yi为xi处的函数估计值n=length(x);m=length(y);%输入的插值点与它的函数值应有相同的个数if n=merror(The lengths of X and Y must be equal!);return;endp=zeros(1,n); %生成n个零元素的行矩阵for k=1:nt=ones(1,n); %生成n个1的行矩阵for j=1:nif j=k % 输入的插值节点必须互异if abs(x(k)-x(j)x=1 2 3 4 5;y=1 8 27 64 125;xi=5;yi=Lagrange(x, y, xi)yi =1252已知0.32 0.34 0.360.314567 0.333487 0.352274求解的近似值用Matlab求解x=0.32 0.34 0.36;y=0.314567 0.333487 0.352274;xi=0.3367; yi=Lagrange(x, y, xi)yi = 0.3304三、复化求积公式3.1算法原理如果积分区间比较大，直接地使用上述求积公式，精度难以保证。在高等教学的学习中，曾经介绍过梯形法求积公式，它的几何意义是将积分区间a, b分成n个小区间，对f (x)用分段线性插值，然后积分。类似地，也可以对f (x)用分段抛物插值。通常采取的办法是：等分求积区间，比如取步长，分a, b为n等分，分点为k = 0, 1, 2, n 在区间 xk, xk+1上使用以上求积公式求得Ik 取和值，作为整个区间上的积分近似值。这种求积方法称为复化求积方法。复化梯型公式3.2实验目的和要求1.熟悉复合梯形求积的基本原理2.通过实例掌握用复合梯形求积的方法3.能编程实现复合梯形求积3.3算法程序function =simpson()min=input(输入区间下限：); %提示输入区间下限max=input(输入区间上限：); %提示输入区间上限n=input(输入等分个数：); %提示输入等分个数syms x;f=input(输入被积函数：); %提示输入被积函数h=(max-min)/n; %求步长b=zeros(1,n);c=zeros(1,n);xa=subs(f,x,min); %求区间下限的函数值m=-xa;u=0;for k=1:n b(k)=min+(k-1)*h; m=m+subs(f,x,b(k);c(k)=b(k)+h/2; u=u+subs(f,x,c(k);endxb=subs(f,x,max); %求区间上限的函数值Sn=(h/6)*(xa+4*u+2*m+xb); %用复化Simpson求积分Sn3.4实例 xk01/81/43/81/25/83/47/81f (xk)43.938463.764703.506853.200002.876402.460002.265492计算积分这个问题有明显的答案取n = 8用复化梯形公式四 runge kutta法4.1算法原理四阶Runger-Kutta法是用, ,和的加权平均值来近似. （4.1）其中，和形如， （4.2） 通过与阶的泰勒级数方法的系数匹配，使得局部误差为，龙格和库塔得出了如下方程组： 必须补充两个条件才能求解，最有用的选择是 （4.3）其余变量的解为： （4.4）将式（4.4）和式（4.5）中的值代入式（4.2）和式（4.1）。得到标准的四阶Runger-Kutta法， 利用 （4.5）生成近似值序列，其中 ， （4.6）综上所述，标准的四阶Runge-Kutta公式为：4.2实验目的与要求1 熟悉runge kutta法的基本原理； 2能编程实现runge kutta法，并求解常微分方程的值；3运用Matlab编程，根据实例中给定的常微分方程，以及初始值，求其在某一点的近似值 4.3算法程序functionx,y=Runge_Kutta4(ydot_fun,x0,y0,h,N)%标准四阶Runge_4方法，其中，%ydot_fun为一阶微分方程的函数；%x0,y0为初始条件；%h为区间步长；%N为区间的个数；%x为Xn构成的向量；%y为Yn构成的向量。x=zeros(1,N+1); %生成N+1个零元素的向量y=zeros(1,N+1); %生成N+1个零元素的向量x(1)=x0;y(1)=y0; %分别把初始值x0，y0赋给向量x,y中的第一个元素for n=1:N %n的取值1，2nx(n+1)=x(n)+h; %x中的元素每次加步长k1=h*feval(ydot_fun,x(n),y(n);k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1);k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(ydot_fun,x(n)+h,y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end4.4实例1求解初始问题的解在处的近似值，取步长。ydot_fun=inline(x+y,x,y);x,y=Runge_Kutta4(ydot_fun,0,1,0.1,1)x = 0 0.1000y = 1.0000 1.11032求解初始问题的解在处的近似值，取步长。ydot_fun=inline(-y,x,y);x,y=Runge_Kutta4(ydot_fun,0,1,0.1,1)x = 0 0.1000y = 1.0000 0.9048五、数值方法实际应用5.1实例地球卫星轨道是一个椭圆,椭圆的周长计算公式为,其中a是椭圆的半长轴,c是地球中心与轨道中心(椭圆中心)的距离,记h为近地点距离,h为远地点距离, R=6371(km)为地球半径,则.我国第一颗人造卫星近地点距离h=439(km),远地点距离h=2384 (km),利用复合Simpson求卫星轨道的周长。5.2程序设计：#include#include#define h 439#define H 2384#define R 6371#define PI 3.1415926void main() double f(double a,double c,double x);int n,i; coutn;double a,c,l;l=PI/(2*n);a=(2*R+H+h)/2;c=(H-h)/2;double s1=0,s2=0,s;for(i=1;in;i+) s1=s1+f(a,c,l*i); for(i=1;i=n;i+) s2=s2+f(a,c,l*(i-0.5);s=4*a*(PI/(12*n)*(f(a,c,0)+f(a,c,PI/2)+2*s1+4*s2);couts=s;double f(double a,double c,double x) return sqrt(1-c*c*sin(x)*sin(x)/(a*a); /*请输入分成的份数 n:20s=48704.5*/心得体会通过两周的数值分析课程设计的学习，我熟练的运用MATLAB软件实现用Gauss消元法求解线性方程组的解,用 Lagrange插值求解在某点值,用二分法求解非线性方程的解，用 runge kutta法求解以阶微分方程在某点的近似值 , 用Jacobi迭代求解线性代数方程组的解.刚开始在编写MATLAB的过程中我遇到了很多问题，开始以为MATLAB的编程语言比较简单很快能上手，但是在实际操作的过程中发现还是有一定的难度，通过前期对C语言的学习，近期对有关知识的理论学习及自学MATLAB的编程语言，总算完成了这份课程设计。虽然选择的课题不是最难的，但是通过这次实习，能促进我进一步学习使用MATLAB在数值分析方面的知识，也算是收获之一。从开始做这篇