偏微分方程实验1

2
1、每种方法的迭代公式：
（1）显式Euler法
p[i]=p[i-1]+h*(*f)(x[i-1],p[i-1]);
（2）隐式梯形法
p[i]=p[i-1]+h*(*f)(x[i-1],p[i-1]);;
（3）二级二阶Runge-Kutta
K1=(*f)(x[i-1],p[i-1]);
K2=(*f)(x[i-1]+h/2,p[i-1]+h*K1/2);
double *p=new double[N];
p[0]=y0;
double t;
for(i=1;i<N;i++)
{
t=p[i-1]+h*fun((i-1)*h,p[i-1]);
p[i]=p[i-1]+0.5*h*(fun((i-1)*h,p[i-1])+fun(i*h,t));
}
return p;
{
return exp(-2*x)-2*x+1;
}
double* Exact(double x0, double h) //精确解
{
int i;
double *p=new double[N];
for(i=0;i<N;i++) p[i]=p_fun(x0+i*h);
return p;
}
double* Euler(double h,double y0) //显式欧拉法
《偏微分方程数值解》
课 程 实 验 报 告（一）
实验名称：
常微分方程初值问题求解
姓 名：
周要
学 号：
指导教师：
王权锋 老师
完成日期：
2019年3月14日
给定初值问题（注：共两次上机时间4学时）
其精确解为 。取h=0.1,分别用显式Euler法、隐式梯形法、二级二阶Runge-Kutta、三级三阶Runge-Kutta、四级四阶Runge-Kutta计算数值解，并与精确解比较。
scanf("%lf",&h);
//调用函数求解
p0=Exact(0,h);
p1=Euler(h,y0);
p2=Y_Euler(h,y0);
p3=Two_R_K(h,y0);
p4=Three_R_K(h,y0);
p5=Four_R_K(h,y0);
//输出
printf("\n x精确解欧拉方法隐式梯形法二阶R_K三阶R_K四阶R_K\n") ;
{
int i;
double *p=new double[N];
p[0]=y0;
for(i=1;i<N;i++) p[i]=p[i-1]+h*fun((i-1)*h,p[i-1]);
return p;
}
double* Y_Euler(double h,double y0) //隐式梯形法
{
int i;
{
k1=fun((i-1)*h,p[i-1]);
k2=fun((i-0.5)*h,p[i-1]+0.5*h*k1);
k3=fun(i*h,p[i-1]+h*(-k1+2*k2));
p[i]=p[i-1]+h/6*(k1+4*k2+k3);
}
return p;
}
double* Four_R_K(double h,double y0) //四阶四级Runge-Kutta
k4=fun(i*h,p[i-1]+h*k3);
p[i]=p[i-1]+h/6*(k1+2*k2+2*k3+k4);
}
return p;
}
//主函数
int main(void)
{
double x,*p0,*p1,*p2,*p3,*p4,*p5;
double h,y0=2;
int i;
//输入步长
printf("请输入步长：");
p[i]=p[i-1]+h*(K1+2*K2+2*K3+K4)/6;
（5）四级四阶Runge-Kutta
K1=(*f)(x[i-1],p[i-1]);
K2=(*f)(x[i-1]+h/2,p[i-1]+h*K1/2);
K3=(*f)(x[i-1]+h/2,p[i-1]+h*K2/2);
K4=(*f)(x[i-1]+h,p[i-1]+h*K3);
学生实验 心得
1、将循环放在主程序中，每次循环调用一次各个算法程序，分别返回p[i]并输出比较方便。
2、高阶龙格库塔公式的精确度更高。
学生（签名）：
年月日
指导
教师
评语成绩评定：指导教（签名）：年 月 日for(i=0;i<N;i++)
printf("%6f %10lf %10lf %10lf %10lf %10lf %10lf\n",h*i,p0[i],p1[i],p2[i],p3[i],p4[i],p5[i]);
}
4
编译运行程序，得到下图
可以看出，欧拉方法的精确度为 ，隐式梯形法和二级二阶Runge-Kutta的精确度为 ，三级三阶Runge-Kutta的精确度为 ，四级四阶Runge-Kutta的精确度为 。
p[i]=p[i-1]+h*k2;
}
return p;
}
double* Three_R_K(double h,double y0) //三阶三级Runge-Kutta
{
int i;
double *p=new double[N];
p[0]=y0;
double k1,k2,k3;
for(i=1;i<N;i++)
p[i]=p[i-1]+h*K2;
（4）三级三阶Runge-Kutta
K1=(*f)(x[i-1],p[i-1]);
K2=(*f)(x[i-1]+h/2,p[i-1]+h*K1/2);
K3=(*f)(x[i-1]+h/2,p[i-1]+h*K2/2);
K4=(*f)(x[i-1]+h,p[i-1]+h*K3);
p[i]=p[i-1]+h*(K1+2*K2+2*K3+K4)/6;
3
#include<stdio.h>
#include<math.h>
#define N 11
double fun(double x,double y) //微分方程
{
return -2*y-4*x;
}
double p_fun(double x) //原函数
}
double* Two_R_K(double h,double y0) //二阶二级Runge-Kutta
{
int i;
double *p=new double[N];
p[0]=y0;
double k1,k2;
for(i=1;i<N;i++)
{
k1=fun((i-1)*h,p[i-1]);
k2=fun((i-0.5)*h,p[i-1]+0.5*h*k1);
{
int i;
double *p=new double[N];
p[0]=y0;
double k1,k2,k3,k4;
for(i=1;i<N;i++)
{
k1=fun((i-1)*h,p[i-1]);
k2=fun((i-0.5)*h,p[i-1]+0.5*h*k1);
k3=fun((i-0.5)*h,p[i-1]+0.5*h*k2);