我尝试编写一个代码来求解第二阶ode，但它以某种方式无法正如我的预期工作。
方程为2y" 5y' 3y = 0，y(0(= 3和y'(0(= -4
最终答案是y(x(= 3-4x＆lt;＆lt;编辑这是错误的
因此y(1(= -1和y'(1(= -4
但是程序输出为y(1(= 0.81414 y'(1(= -1.03727＆lt;＆lt;正确
请帮助！
谢谢！
#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <math.h>
double ddx(double x,double y,double z)
{
double dx = (-5*z-3*y)/2;
return dx;
}
double test2ndorder(double x0, double y0, double z0, double xt, double h)
{
int n = (int)((xt - x0) / h);
double k1, k2, k3, k4;
double l1, l2, l3, l4;
double x = x0;
double y = y0;
double z = z0;
for (int i = 1; i <= n; i++)
{
k1 = h * z;
l1 = h * ddx(x, y, z);
k2 = h * (z + 0.5*l1);
l2 = h * ddx(x + 0.5*h, y + 0.5*k1, z + 0.5*l1);
k3 = h * (z + 0.5*l2);
l3 = h * ddx(x + 0.5*h, y + 0.5*k2, z + 0.5*l2);
k4 = h * (z + l3);
l4 = h * ddx(x + h, y + k3, z + l3);
y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4);
z = z + (1.0 / 6.0)*(l1 + 2 * l2 + 2 * l3 + l4);
x = x + h;
std::cout << y << " ";
std::cout << z << "n";
}
return y;
}
int main()
{
double x0, y0, z0, x, y, z,h;
x0 = 0;
x = 1;
y0 = 3;
z0 = -4;
h =0.01;
y = test2ndorder(x0, y0, z0, x, h);
std::cout << y;
}
在 test2ndorder中的for循环中，您使用常量的 x0， y0和 z0值，因此每次迭代中计算的值都不会更改。您还计算了y的新值，并且永远不会使用它。
看来您应该使用x，y和z。
2y" + 5y' +3y = 0 ,y(0) = 3 and y'(0) = -4
具有特征根，是
的解决方案(2r)^2 + 5*(2r) + 6 = 0 <==> r = -1 or r = -1.5
使解决方案为
y(x) = A*exp(-x)+B*exp(-1.5*x)
3 = y(0) = A + B
-4 = y'(0) = -A - 1.5*B
，可以给出
的B = 2和A = 1 y(1) = 0.814139761468302
y'(1) = -1.03726992161673
这是关于您得到的。
C 语言中的代码。在此实现中，使用了以下初始条件，x0 = 0，y(x0(= 3和u(x0(= -4。
// Autor : Carlos Eduardo da Silva Lima
// Tema : Runge-Kutta de quarta ordem
// Linguagem: C++ (ANSI)
// IDE : Online GBD
// Data : 19/11/2022
#include <iostream>
#include <cmath>
#define N 10000
float edo1(float x, float y, float u);
float edo2(float x, float y, float u);
void RK4(float x0, float y0, float u0, float h);
int main()
{
RK4(0, 3, -4, 1E-3);
return 0;
}
float edo1(float x, float y, float u){
return u;
}
float edo2(float x, float y, float u){
return ((-3*y-5*u)/2);
}
void RK4(float x0, float y0, float u0, float h){
float x[N], y[N], u[N];
x[0] = x0;
y[0] = y0;
u[0] = u0;
int i = 0;
do{
float k11 = h*edo1(x[i],y[i],u[i]);
float k12 = h*edo2(x[i],y[i],u[i]);
float k21 = h*edo1(x[i]+(h/2),y[i]+(k11/2),u[i]+(k12/2));
float k22 = h*edo2(x[i]+(h/2),y[i]+(k11/2),u[i]+(k12/2));
float k31 = h*edo1(x[i]+(h/2),y[i]+(k21/2),u[i]+(k22/2));
float k32 = h*edo2(x[i]+(h/2),y[i]+(k21/2),u[i]+(k22/2));
float k41 = h*edo1(x[i]+h,y[i]+k31,u[i]+k32);
float k42 = h*edo2(x[i]+h,y[i]+k31,u[i]+k32);
y[i+1] = y[i] + ((k11+2*(k21+k31)+k41)/6);
u[i+1] = u[i] + ((k12+2*(k22+k32)+k42)/6);
x[i+1] = x[i] + h;
i++;
}while(i<=(N-1));
// Saída
int j = 0;
do{
std::cout << "x = " << x[j] << " | y = " << y[j] << " | u = " << u[j] << std::endl;
j++;
}while(j<=(N-1));
}
我也用Python语言进行了实现。在这里，我在scipy， scipy.integrate.odeint 中使用了一种集成方法，用于EDO，并在Scipy中开发和记录。如果您不同意，或者发现任何逻辑错误，请告诉我。直到：(。
# Autor : Carlos Eduardo da Silva Lima
# Tema : Runge-Kutta de quarta ordem
# Linguagem: Python
# IDE : Google Colab
# Data : 19/11/2022
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# Condições iniciais
x_inicial = 0.0
x_final = 10.0
y0 = 3
u0 = -4
N = 10000
# Equações diferenciais na primeira ordem
def edo(r,x):
y = r[0]
u = r[1]
return np.array([u, (-3*y-5*u)/2])
# Aplicação do métodos para resolução das edo´s declaradas acima emonjunto com as condições
x = np.linspace(x_inicial, x_final, N)
r0 = np.array([y0,u0])
# Resolução
sol = odeint(edo, r0, x)
# Grvandos os resultados em novas variáveis
y = sol[:,0]
u = sol[:,1]
# Plot
plt.style.use('dark_background')
plt.figure(figsize=(7,7))
plt.plot(x,y,'m-', linewidth = 3.5)
plt.grid(color='y', linestyle='-', linewidth=0.05)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Odeint python")
plt.show()