c - C 上的 Runge-Kutta 积分

Question

所以我的代码有效，只是出于某种原因，我的代码总是运行两个 if 语句(两个 y 方程，无论我为第一个 fprintf 问题输入哪个数字)。此外，t,y 列总是比 t,y2 列长得多(编辑，即如果我输入 funt=1,t0=0,t1=10,inity=1,dt=0.1，我的 t,y 列有 100行和我的 t,yexact 列有 10 行)。此外，t,y 部分的 t 列对于第一个 if 语句(y 方程)全为零，对于第二个 if 语句(y 方程)全为 1。任何帮助将不胜感激。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double rk4(double(*f)(double,double),double,double,double);
double user_fun_1(double,double);
double user_fun_2(double,double);

int main (void) {
  int funct;
  double *y,t,y2;
  double t0,t1,dt,inity,C;
  int i,n;
  fprintf(stdout, "Which function do you wish to integrate?\n(1)f(t,y)=t\n(2)f(t,y)=-y\n");
  fscanf(stdin,"%d",&funct);
  fprintf(stdout,"What is the initial time?\n");
  fscanf(stdin,"%lf",&t0);
  fprintf(stdout,"What is the final time?\n");
  fscanf(stdin,"%lf",&t1);
  fprintf(stdout,"What is the initial value of y?\n");
  fscanf(stdin,"%lf",&inity);
  fprintf(stdout,"what is the step size?\n");
  fscanf(stdin,"%lf",&dt);
  if (dt>(t1-t0)){
    printf("The stepsize has to be less than or equal to the total time\n");
    return EXIT_FAILURE;}
  n = 1 + (t1-t0)/dt;
  y = (double *)malloc(n*sizeof(double));
  if (funct=1){
    C=inity-(t0*t0*0.5);
    printf("%s%25s\n","t","y");
    for (y[0]=inity,i=1;i<n;i++){
      y[i]=rk4(user_fun_1,dt,t0 + dt * (i-1),y[i-1]);
      printf("%lf%25lf\n",t,y[i]);}
    printf("%s%25s\n","t","y exact");
    for(i=0;i<n;i+=t1){
      t=t0+dt*i;
      y2=(t*t*0.5) + C;
      printf("%lf%25lf\n",t,y2);}}

    if (funct=2){
      C=inity/(exp(-t0));
      printf("%s%25s\n","t","y");
      for(y[0]=inity,i=1;i<n;i++){
        y[i]=rk4(user_fun_2,dt,t0+dt*(i-1),y[i-1]);
        printf("%lf%25lf\n",t,y[i]);}
        printf("%s%25s\n","t","y exact");
      for(i=t0;i<n;i+=t1){
        t=t0+dt*i;
        y2=C*exp(-t);
        printf("%lf%25lf\n",t,y2);}}

  free(y);
  return EXIT_SUCCESS;}

double user_fun_1(double t,double y){
  return t;}
double user_fun_2(double t, double y){
  return -y;}
double rk4(double(*f)(double, double), double dt,double t,double y){
  double k1,k2,k3,k4;
  k1=dt*f(t,y);
  k2=dt*f(t + dt/2,y + k1/2);
  k3=dt*f(t + dt/2,y + k2/2);
  k4=dt*f(t + dt,y + k3);
  return y + (k1 + 2*k2 + 2*k3 +k4)/6;}

Answer 1

您的if 语句是错误的，它是赋值而不是比较。应该是

if (funct==2)

不

if (funct=2)

第一个 if 语句也是如此。

一个安全的 c 习惯是在 if 语句中交换常量和变量的位置。 if(2==funtc) 有效，而 if(2=funct) 会给出编译时错误。下次尝试打开编译器警告，它可以帮助您发现此类错误。如果您使用的是 gcc，请添加 -Wall -Wextra 标志。

你的 t, y 比 t, y exact 长的原因是你的 n =101，而 t1 = 10 。所以在第一个 for 循环中，你写了 i++，所以它循环了 100 次，但是在第二个 for 循环中，你使用了 i+=t1，所以只循环了11次。

t 始终为 0 的原因是因为您从未为其分配任何初始值，也没有在第一个 for 循环中分配值。