| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 1806 人关注过本帖
标题:偏微分方程的数值模拟的程序,输出的解为0,求修改
只看楼主 加入收藏
busong
Rank: 1
等 级:新手上路
帖 子:3
专家分:0
注 册:2020-6-23
结帖率:0
收藏
已结贴  问题点数:20 回复次数:3 
偏微分方程的数值模拟的程序,输出的解为0,求修改
#include <cstdlib>
#include <cmath>
#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>


using namespace std;

const int Q = 4;     
const int NX =400;  

int e[Q] = {1,-1,2,-2};
double s[Q] ={0.7,19.0/30.0,-11.0/60.0,-3.0/20.0};
double u[NX+3],u1[NX+3],u2[NX+3],Du[NX+3],u0[NX+3],ut[NX+3],utt[NX+3],ut0[NX+3],uxx[NX+3],f[NX+3][Q],F[NX+3][Q];
double g[NX+3],sour_xx[NX+3];
double c,dx,dt,tau,error,maxu;
double x_start,x_end;
double a;

int i,j,ip;

void init();
double feq(int k,double u,double ut);
void evolution(double t);
void output(int m);
void Error();
double exact_solution(double x,double t);
double exact_solution_ut(double x,double t);
inline double fai_one(double u);
inline double sour(double u);
inline double sour_t(double u,double ut);
inline double sour_tt(double u,double ut,double utt);
void u_time(double t);
void prepare();
void max();


int main()
{
    int m=1;
    int n=1;
    double t=0;
    x_start=-2.0;
    x_end=2.0;
    dx=(x_end-x_start)/NX;
    c=1.0;
    dt=dx/c;

    ostringstream name;
    name<<"GRE for different tau"<<".dat";
    ofstream out(name.str().c_str());
    out<<"dx="<<dx<<",dt="<<dt<<",c="<<c<<",a="<<a<<endl;
    out<<"s[Q]="<<s[0]<<","<<s[1]<<","<<s[2]<<","<<s[3]<<endl;

    tau=0.5;
    for(m=1;m<=20;m++)
    {
        tau=tau+0.001;
        t=0;
        init();
        for(n=1;n<=2000;n++)
        {
            t=n*dt;
            evolution(t);
            if(n%200==0)
            {
                u_time(t);
                max();
                Error();
                out<<tau<<"t"<<n<<"t"
                <<resetiosflags(ios::fixed)<<setiosflags(ios::scientific)<<setprecision(6)
                <<error<<"t"<<maxu<<resetiosflags(ios::scientific)<<endl;
            }
        }
    }
    return 0;
}

void init()
{
    a=1.0;
    double x=0;
    double t0=0;
   
    u_time(t0);
    for(i=0;i<=NX+2;i++)
    {
        x=x_start+double(i-1)*dx;
        u[i]=exact_solution(x,t0);
        u1[i]=u[i];
        u2[i]=u[i];
        ut0[i]=4.0/cosh(x);
        utt[i]=0.0;
        
        for(j=0;j<Q;j++)
        {
            f[i][j]=feq(j,u[i],ut[i]);
            F[i][j]=0.0;
        }
        
    }
    output(0);
}

void prepare()
{
    for(i=0;i<=NX+2;i++)
    {
        g[i]=-a*a*u[i]*sin(u[i]);
    }

    for(i=2;i<=NX;i++)
    {
        sour_xx[i]=(-g[i+2]+16.0*g[i+1]-30.0*g[i]+16.0*g[i-1]-g[i-2])/(12.0*dx*dx);
    }

    sour_xx[0]=(-35.0*g[4]+104.0*g[3]-114.0*g[2]+56.0*g[1]-11.0*g[0])/(12.0*dx*dx);
    sour_xx[1]=(-35.0*g[5]+104.0*g[4]-114.0*g[3]+56.0*g[2]-11.0*g[1])/(12.0*dx*dx);

    sour_xx[NX+1]=(35.0*g[NX+1]-104.0*g[NX]+114.0*g[NX-1]-56.0*g[NX-2]+11.0*g[NX-3])/(12.0*dx*dx);
    sour_xx[NX+2]=(35.0*g[NX+2]-104.0*g[NX+1]+114.0*g[NX]-56.0*g[NX-1]+11.0*g[NX-2])/(12.0*dx*dx);
   
    for(i=0;i<=NX+2;i++)
    {
        ut[i]=(u[i]-u1[i])/dt;
        utt[i]=(u[i]+u2[i]-2.0*u1[i])/(dt*dt);
    }
}

void evolution(double t)
{   
    prepare();
    int i=0;
    int j=0;
    int ip=0;

    for(i=2;i<NX+1;i++)
    {
        for(j=0;j<Q;j++)
        {
            double f_s1,f_s2,f_s3,f_s4;
            ip=i-e[j];
            f_s1=dt*sour(u[ip])*s[j];
            f_s2=dt*dt*sour_t(u[ip],ut[ip])*0.5*s[j];
            f_s3=dt*dt*dt*sour_tt(u[ip],ut[ip],utt[ip])*s[j]/6.0;
            f_s4=dt*dt*(2.0*tau*tau-2.0*tau+0.25)*sour_xx[ip]*s[j]/(tau-0.5);

            F[i][j]=f[ip][j]+(feq(j,u[ip],ut[ip])-f[ip][j])/tau+f_s1+f_s2+f_s3+f_s4;
        }
    }

    for(i=0;i<=NX+2;i++)
    {
        u2[i]=u1[i];
        u1[i]=u[i];
    }

    for(i=2;i<NX+1;i++)
    {
        ut[i]=0.0;
        for(j=0;j<Q;j++)
        {
            f[i][j]=F[i][j];
            ut[i]+=f[i][j];
        }
        
    }
    ut[1]=exact_solution_ut(x_start,t);
    ut[0]=exact_solution_ut(x_start-dx,t);
    ut[NX+1]=exact_solution_ut(x_end,t);
    ut[NX+2]=exact_solution_ut(x_end+dx,t);

    f[1][0]=feq(0,u[1],ut[1])+(f[2][0]-feq(0,u[2],ut[2]));
    f[1][2]=feq(2,u[1],ut[1])+(f[2][2]-feq(2,u[2],ut[2]));
    f[0][2]=feq(2,u[0],ut[0])+(f[1][2]-feq(2,u[1],ut[1]));
    f[NX+1][1]=feq(1,u[NX+1],u[NX+1])+(f[NX][1]-feq(1,u[NX],ut[NX]));
    f[NX+1][3]=feq(3,u[NX+1],u[NX+1])+(f[NX][3]-feq(3,u[NX],ut[NX]));
    f[NX+2][3]=feq(3,u[NX+2],u[NX+2])+(f[NX+1][3]-feq(3,u[NX+1],ut[NX+1]));
}

double feq(int k,double u,double ut)
{
    double feq=0;
    switch(k)
    {
    case 0:feq=(4.0*ut-fai_one(u))/6.0;
        return feq;
    case 1:feq=(4.0*ut-fai_one(u))/6.0;
        return feq;
    case 2:feq=(fai_one(u)-ut)/6.0;
        return feq;
    case 3:feq=(fai_one(u)-ut)/6.0;
        return feq;
    default:
        return 0;
    }
}

void output(int m)
{
    double x=0;
    int i=0;

    ostringstream name;
    name<< "D1Q4_"<<m<<".dat";
    ofstream out(name.str().c_str());
    out << "Title = "heat conduction ep"n"
        << "VARIABLES = "X","U","U0"n"
        << "tau="<<tau<<" "<<NX+1<<"POINTs"<<" "<<"dx="<<dx<<" "<<"dt="<<dt<<" "<<endl;

    for(i=0;i<=NX+2;i++)
    {
        x=x_start+double(i-1)*dx;
        out << x << "t" <<u[i]<<"t"<<u0[i]<<endl;
    }
}

void Error()
{
    double temp1,temp2;
    temp1=0;
    temp2=0;
   
    for(i=1;i<=NX+1;i++)
    {
        temp1+=fabs(u[i]-u0[i]);
        temp2+=fabs(u0[i]);
    }
        error=temp1/temp2;
}

double exact_solution(double x,double t)
{
    double u;
    u=4.0*atan(t/cosh(x));
    return u;
}

 double exact_solution_ut(double x,double t)
{
    double u;
    u=4.0/(cosh(x)+t*t/cosh(x));
    return u;
}

inline double fai_one(double u)
{
    double fai_one;
    fai_one=u*pow(a,2.0)/(dt*(tau-0.5)*c*c);
    return fai_one;
}

inline double sour(double u)
{
    double sour;
    sour=-sin(u);
    return sour;
}

inline double sour_t(double u,double ut)
{
    double sour_t;
    sour_t=-ut*cos(u);
    return sour_t;
}

inline double sour_tt(double u,double ut,double utt)
{
    double sour_tt;
    sour_tt=-utt*cos(u)-ut*sin(u)*ut;
    return sour_tt;
}

void u_time(double t)
{
    double x=0;

    for(i=0;i<=NX+2;i++)
    {
        x=x_start+double(i-1)*dx;
        u0[i]=exact_solution(x,t);
    }
}

void max()
{
    double abu;
    maxu=0;
   
    for(i = 0;i <= NX+2 ;i++)
    {
        abu=fabs(u[i]-u0[i]);
        if (abu>maxu)
            maxu=abu;
    }
}
输出的解为0,请问怎么修改
搜索更多相关主题的帖子: void for ip int double 
2020-06-23 14:05
busong
Rank: 1
等 级:新手上路
帖 子:3
专家分:0
注 册:2020-6-23
收藏
得分:0 
有没有人会呀
2020-06-24 09:43
ditg
Rank: 10Rank: 10Rank: 10
等 级:贵宾
威 望:16
帖 子:852
专家分:1937
注 册:2014-4-10
收藏
得分:20 
学妹啊,俺是真的忘了,就算不会也得知会一声……

换人,下一位!

梦想拥有一台龙芯3A-4000
2020-06-24 10:10
busong
Rank: 1
等 级:新手上路
帖 子:3
专家分:0
注 册:2020-6-23
收藏
得分:0 
怎么能删掉我发的这个帖子
2021-03-02 21:57
快速回复:偏微分方程的数值模拟的程序,输出的解为0,求修改
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.017312 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved