第四周3.中的三个根。

，在求8] [0,041.76938.7911.1-)(23=-+=x x x x ffunction y=mj()for x0=0:0.01:8x1=x0^3-11.1*x0^2+38.79*x0-41.769;if (abs(x1)<1.0e-8)x0endend4.分别用简单迭代法、埃特金法、牛顿法求解方程，并比较收敛性与收敛速度（ε分别取10-3、10-5、10-8）。

简单迭代法：function y=jddd(x0)x1=(20+10*x0-2*x0^2-x0^3)/20;k=1;while (abs(x1-x0)>=1.0e-3)x0=x1;x1=(20+10*x0-2*x0^2-x0^3)/20;k=k+1;endx1k埃特金法：function y=etj(x0)x1=(20-2*x0^2-x0^3)/10;x2=(20-2*x1^2-x1^3)/10;x3=x2-(x2-x1)^2/(x2-2*x1+x0);k=1;while (abs(x3-x0)>=1.0e-3)x0=x3;x1=(20-2*x0^2-x0^3)/10;x2=(20-2*x1^2-x1^3)/10;x3=x2-(x2-x1)^2/(x2-2*x1+x0);k=k+1;end2,020102)(023==-++=x x x x x fx3k牛顿法：function y=newton(x0)x1=x0-fc(x0)/df(x0);k=1;while (abs(x1-x0)>=1.0e-3)x0=x1;x1=x0-fc(x0)/df(x0);k=k+1;endx1kfunction y=fc(x)y=x^3+2*x^2+10*x-20;function y=df(x)y=3*x^2+4*x+10;第六周1.解例6-4（p77)的方程组，分别采用消去法（矩阵分解）、Jacobi迭代法、Seidel迭代法、松弛法求解，并比较收敛速度。

消去法：x=a\d或[L,U]=lu(a);x=inv(U)inv(L)dJacobi迭代法：function s=jacobi(a,d,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);C=inv(D);B=C*(L+U);G=C*d;s=B*x0+G;n=1;while norm(s-x0)>=1.0e-8x0=s;s=B*x0+G;n=n+1;endnSeidel迭代法：function s=seidel(a,d,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);C=inv(D-L);B=C*U;G=C*d;s=B*x0+G;n=1;while norm(s-x0)>=1.0e-5x0=s;s=B*x0+G;n=n+1;endn松弛法：function s=loose(a,d,x0,w)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);C=inv(D-w*L);B=C*((1-w)*D+w*U);G=w*C*d;s=B*x0+G;n=1;while norm(s-x0)>=1.0e-8x0=s;s=B*x0+G;n=n+1;endn2.练习MATLAB的常用矩阵语句，就龙格现象函数（p88）练习插值语句interp, spline，并比较。

3.测得血液中某药物浓度随时间的变化值为:求t=0.45, 1.75, 5.0, 6.0 时的浓度C.分别用n=4,5,9的拉格朗日插值计算；并用样条函数插值计算,并比较结果。

拉格朗日插值:function s=lagr(n)x=[0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0];y=[19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88];x0=[0.45 1.75 5.0 6.0];m=length(x0);for i=1:mD=abs(x-x0(i));I=1;while I<=n+1for a=1:length(x)if D(a)==min(D)c(I)=a;D(a)=max(D)+1;breakendendI=I+1;endb=sort(c);z=x0(i);t=0.0;for k=1:length(b)u=1.0;for j=1:length(b)if j~=ku=u*(z-x(b(j)))/(x(b(k))-x(b(j)));endendt=t+u*y(b(k));ends(i)=t;end样条函数差值：Interp1(x,y,x0,’spline’)Spline(x,y,x0)第八周1.给定某药物浓度随时间的变化值(作业3)，1）分别采用样条函数和三点公式（设h=0.1)求结点处的导数值，并比较结果。

2）求该时间段的平均浓度（定步长S法）样条函数：x=[0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0];y=[19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88];pp=csape(x,y,'not-a-knot');df=fnder(pp);df1=ppval(df,x)三点公式：function df=sandian()t=[0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0];c=[19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88];h=0.1;n=length(t);for i=1:nx0=t(i);y0=c(i);y1=spline(t,c,x0+h);y2=spline(t,c,x0+2*h);y3=spline(t,c,x0-h);y4=spline(t,c,x0-2*h);switch icase 1df(i)=(-3*y0+4*y1-y2)/(2*h);case ndf(i)=(y4-4*y3+3*y0)/(2*h);otherwisedf(i)=(y1-y3)/(2*h);endendend平均浓度：function averagec=simpson()t=[0.25 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0];c=[19.30 18.15 15.36 14.10 12.89 9.32 7.55 5.24 3.86 2.88]; m=(t(1)+t(10))/2;y=spline(t,c,m);averagec=(c(1)+4*y+c(10))/6;end2.练习MATLAB常用的trapz, quad, quadl等语句。

计算：x=0:8;y=1./(sqrt(2.*pi)).*exp(-(x-4).^2./2);z=trapz(x,y)function y=jifen(x)y=1./(sqrt(2.*pi)).*exp(-(x-4).^2./2);q1=quad('jifen',0,8,1.0e-8)q2=quadl('jifen',0,8,1.0e-8)3.采用变步长经典R-K法, ode23, ode45计算例9-5，并作比较。

变步长经典R-K法：（可能有问题）function z=jdrk(m)x0=[25 2]';a=0;b=15;n=length(x0);z=zeros(n,m);k1=zeros(n,1);k2=zeros(n,1);k3=zeros(n,1);k4=zeros(n,1);t=a;x=x0;x2=zeros(n,1);x3=x2;x4=x2;h=choose(m);m1=15/h+1;for k=1:m1k1=prey(t,x);for i=1:nx2(i)=x(i)+1/2*h*k1(i);endk2=prey(t+h/2,x2);for i=1:nx3(i)=x(i)+1/2*h*k2(i);endk3=prey(t+h/2,x3);for i=1:nx4(i)=x(i)+h*k3(i);endk4=prey(t+h,x4);for i=1:nx(i)=x(i)+h/6*(k1(i)+2*k2(i)+2*k3(i)+k4(i)); z(i,k)=x(i);endt=t+h;endh1=length(z);t2=[a:(b-a)/(h1-1):b];plot(t2,z)gtext('x1(t)')gtext('x2(t)')function h=choose(n)h=15/(n-1);t0=0;x0=[25 2]';k11=prey(t0,x0);k21=prey(t0+h/2,x0+h/2*k11);k31=prey(t0+h/2,x0+h/2*k21);k41=prey(t0+h,x0+h*k31);x1=x0+h/6*(k11+2*k21+2*k31+k41);k12=prey(t0,x0);k22=prey(t0+h/4,x0+h/4*k12);k32=prey(t0+h/4,x0+h/4*k22);k42=prey(t0+h/2,x0+h/2*k32);x2=x0+h/12*(k12+2*k22+2*k32+k42);if abs(x2-x1)<1.0e-5while abs(x2-x1)<1.0e-5h=h*2;k11=prey(t0,x0);k21=prey(t0+h/2,x0+h/2*k11);k31=prey(t0+h/2,x0+h/2*k21);k41=prey(t0+h,x0+h*k31);x1=x0+h/6*(k11+2*k21+2*k31+k41);k12=prey(t0,x0);k22=prey(t0+h/4,x0+h/4*k12);k32=prey(t0+h/4,x0+h/4*k22);k42=prey(t0+h/2,x0+h/2*k32);x2=x0+h/12*(k12+2*k22+2*k32+k42);endh=h/2;elsewhile abs(x2-x1)>=1.0e-5h=h/2;k11=prey(t0,x0);k21=prey(t0+h/2,x0+h/2*k11);k31=prey(t0+h/2,x0+h/2*k21);k41=prey(t0+h,x0+h*k31);x1=x0+h/6*(k11+2*k21+2*k31+k41);k12=prey(t0,x0);k22=prey(t0+h/4,x0+h/4*k12);k32=prey(t0+h/4,x0+h/4*k22);k42=prey(t0+h/2,x0+h/2*k32);x2=x0+h/12*(k12+2*k22+2*k32+k42);endendfunction xdot=prey(t,x)r=1;a=0.1;b=0.5;c=0.02;xdot=[r-a*x(2) 0;0 -b+c*x(1)]*x;ode23, ode45：[t,x]=ode23('prey',[0:0.1:15],[25 2]);plot(t,x)[t,x]gtext('x1(t)')gtext('x2(t)')[t,x]=ode45('prey',[0:0.1:15],[25 2]);plot(t,x)[t,x]gtext('x1(t)')gtext('x2(t)')第十周1．熟悉常用的概率分布、概率密度函数图、分位点。