matlab程序设计实践

一、 1、高斯求积法

（1）算法说明

为了使求积公式得到较高的代数精度，可以使用到Guass公式：

?1?1f(x)dx??Akf(xk)

k?0n式中Ak为Gauss-Legendre系数。

它可使计算达到2n-1次的代数精度。对任意的求积区间[a,b]，通过变换x=(a-b)x/2+(a+b)/2，可以转化到区间[-1,1]上。此时

?baf(x)dx?b?a1b?aa?bf(t?)dt 2??122（2）流程图

计算步长 h=(b-a)/n s=0 m=0 是 m>n/2-1 否 s=s+h*(guassf(a+h*((1-1/sqrt(3))+2*m))+guassf(a+h*((1+1/sqrt(3))+2*m))) m=m+1 输出s

（3）源程序代码

function s=guass(a,b,n) h=(b-a)/n; s=0.0;

for m=0:(1*n/2-1)

s=s+h*(guassf(a+h*((1-1/sqrt(3))+2*m))+guassf(a+h*((1+1/sqrt(3))+2*m))); end

2、算法应用 求下面的积分值

G??cos(x)dx

01（1）流程图

由高斯余项确定n值,计算得到n=6 建立目标函数guassf y=cos(x) 调用Guass求积函数guass计算目标函数在[0,1]的积分 得出结果

（2）源程序代码 function y=guassf(x) y=cos(x);

guass(0,1,6) ans = 0.8415

二、已知Appolo卫星的运动轨迹（x,y）满足下面方程：

d2xdy?(x??)?(x??)， ?2?x??233dtdtr1r2

d2ydx?y?y??2?y??3 23dtdtr1r2其中，??1/82.45,??1??,r1?(x??)2?y2,r2?(x??)2?y2 ，试在

初值x(0)=1.2,x’(0)=0, y(0)=0, y’(0)=-1.04935371下求解，并绘制Appolo卫星轨迹图。

（1）流程图

建立目标函数appollo(t,x) 给常量mu,lamda及变量r1,r2赋值 令x=[x ；x’； y； y’] 则dx=[x’； x’’； y’ ；y’’] 设定初值 x(0)=1.2，x’(0)=0 y(0)=0,y’(0)=-1.04935371 设置积分的相对误差1e-8 积分限定为 [0，20] 调用常微分方程函数ode45求出数值解 调用绘图函数plot绘出x对y的图形

（2）源程序代码

function dx=appollo(t,x) mu=1/82.45; lamda=1-mu;

r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)+lamda)^2+x(3)^2); dx=[x(2);

2*x(4)+x(1)-lamda*(x(1)+mu)/r1^3-mu*(x(1)-lamda)/r2^3; x(4);

-2*x(2)+x(3)-lamda*x(3)/r1^3-mu*x(3)/r2^3];

x0=[1.2;0;0;-1.04935371]; options=odeset('reltol',1e-8);

[t,y]=ode45('appollo',[0,20],x0,options); plot(y(:,1),y(:,3))

title('Appollo卫星运动轨迹') xlabel('X') ylabel('Y')

（3）运行结果

三、已知：

s?1?2?22?23?...?263分别用循环结构和调用

MATLAB的sum函数求s的值。

（1）流程图