内容发布更新时间 : 2024/11/19 14:46:53星期一 下面是文章的全部内容请认真阅读。

}

break;

case 3: //以下是三次样条函数插值

d[0]=6/h*((f(xi[1])-f(xi[0]))/h-f1(xi[0])); printf(\ for(j=1;j<=N-1;j++)

{ d[j]=6*chashan(xi[j-1],xi[j],xi[j+1]); printf(\ }

d[N]=6/h*(f1(xi[N])-(f(xi[N])-f(xi[N-1]))/h); printf(\

printf(\请用Gauss列主元法解十一阶三对角方程求出M0至M10 !\

printf(\然后把M0至M10输入! 准备好后按1:\scanf(\ if(cchar==1)

{ for(j=0;j<=10;j++)

{ printf(\请输入M%d:\ scanf(\ } }

for(j=0;j<=9;j++) for(i=0;i<=N;i++)

{if(xo1[j]>xi[i] && xo1[j]

help2=S3(xi[i],xi[i+1],M[i],M[i+1],xo1[j]); helpe=help1-help2;

helpb=helpe/help1*100.0;

printf(\ printf(\ printf(\printf(\相对误差=%9.6f%% \ }

if(xo2[j]>xi[i] && xo2[j]

help2=S3(xi[i],xi[i+1],M[i],M[i+1],xo2[j]); printf(\ printf(\ } }

11

附2：综合实验参考 （1） 设计思路：

根据欧拉方法、改进的欧拉方法和经典的4级4阶Runge-Kutta法来求解。 （2） 函数文件： 公用的三个函数文件： fun.m函数文件： function y=fun(x1,x2) if -0.8*x1+0.3*x1*x2~=0

y=(1.2*x2-0.6*x1*x2)/(-0.8*x1+0.3*x1*x2); end

fun1.m函数文件： function y=fun1(x) y=-0.8*x(1)+0.3*x(1)*x(2); fun2.m函数文件： function y=fun2(x) y=1.2*x(2)-0.6*x(1)*x(2); 第一问的三种方法求解文件： euler方法： function euler(n,h) x(1)=2; x(2)=1;

12

t(1)=0.0000; p(1)=x(1); q(1)=x(2); for i=1:n

f(1)=fun1(x); f(2)=fun2(x); y=x+h*f; x=y; t(i+1)=i*h; p(i+1)=x(1); q(i+1)=x(2); end [t',p',q'];

plot(t,p,'r',t,q,'k'); xlabel('时刻t'); ylabel('数量x'); text(0,2,'蓝鲸x1'); text(0,1,'磷虾x2'); 改进的euler方法： function impeuler(n,h) x(1)=2; x(2)=1;

13

t(1)=0.0000; p(1)=x(1); q(1)=x(2); for i=1:n

f(1)=fun1(x); f(2)=fun2(x); g(1)=fun1(x+h*f); g(2)=fun2(x+h*f); x=x+1/2*h*(f+g); t(i+1)=i*h; p(i+1)=x(1); q(i+1)=x(2); end [t',p',q'];

plot(t,p,'r',t,q,'k'); xlabel('时刻t'); ylabel('数量x'); text(0,2,'蓝鲸x1'); text(0,1,'磷虾x2');

经典的4级4阶Runge-Kutta法： function rk4(n,h) x(1)=2;

14

x(2)=1; t(1)=0.0000; p(1)=x(1); q(1)=x(2); for i=1:n

s1(1)=h*fun1(x); s1(2)=h*fun2(x); s2(1)=h*fun1(x+1/2*s1); s2(2)=h*fun2(x+1/2*s1); s3(1)=h*fun1(x+1/2*s2); s3(2)=h*fun2(x+1/2*s2); s4(1)=h*fun1(x+s3); s4(2)=h*fun2(x+s3); x=x+1/6*(s1+2*s2+2*s3+s4); t(i+1)=i*h; p(i+1)=x(1); q(i+1)=x(2); end [t',p',q'];

plot(t,p,'r',t,q,'k'); xlabel('时刻t'); ylabel('数量x');

15