问题描述
我正在尝试实现欧拉方法来求解具有初始条件 y' = x^3 + y^2 的 0 到 2 之间的微分方程 y(0) = 0.5。
首先我设置了 h = 0.1 并且没问题。 y 收敛到 643，....
其次，我将 h = 0.01 和 y 发散到无穷大。
谁能解释为什么会这样？
我的 MATLAB 代码如下：
clc;
clear;
close all;
% dy/dt = f(t,y)
f = @(t,y) t.^3 + y.^2;
a = 0;
ya = 0.5;
b = 2;
h = 0.1;
% h = 0.01;
M = (b-a)/h;
E = euler(f,a,b,ya,M);
function [E] = euler(f,M)
% INPUT - f is the function
% - a and b are the left and right endpoints
% - ya is the initial condition y(a)
% - M is the number of steps
% OUTPUT - E = [T Y] where T is the vector of abscissas and
% Y is the vector of ordinates
h = (b-a) / M;
T = (a:h:b)';
Y = zeros(M+1,1);
Y(1) = ya;
for j = 1:M
Y(j+1) = Y(j) + h*f(T(j),Y(j));
end
E = [T Y];
end
解决方法
发生这种情况是因为精确解在那个区间有一个奇点，随着步长的减小，你越来越接近精确解及其散度。
您可以使用代换 class Author(models.Model):
name = models.CharField(max_length=255)
created_at = models.DateTimeField(auto_now_add = True)
updated_at = models.DateTimeField(auto_now = True)
class Book(models.Model):
title = models.CharField(max_length=255)
contributor = models.ForeignKey(User,related_name = "book_added",on_delete = models.CASCADE)
author = models.ForeignKey(Author,related_name = "book_written",on_delete = models.CASCADE)
created_at = models.DateTimeField(auto_now_add = True)
updated_at = models.DateTimeField(auto_now = True)
class Review(models.Model):
content = models.CharField(max_length=255)
rating = models.IntegerField()
poster = models.ForeignKey(User,related_name="user_reviews",on_delete=models.CASCADE)
book = models.ForeignKey(Book,related_name="book_reviews",on_delete=models.CASCADE)
created_at = models.DateTimeField(auto_now_add = True)
updated_at = models.DateTimeField(auto_now = True)
来更安全地求解这个 Riccati 方程以获得二阶方程 y=-u'/u。 u''(x)+x^3*u(x)=0 的根是 u 的奇点。代码
y给出一个情节
在大约 f = lambda x,u: np.array([u[1],-x**3*u[0]])
u0 = np.array([1,-0.5])
x_ref = np.linspace(0,2,201)
u_ref = odeint(f,u0,x_ref,tfirst=True,atol=1e-8,rtol=1e-10)
x_euler=np.linspace(0,21)
dx = x_euler[1]-x_euler[0]
u_euler = np.array(len(x_euler)*[u0])
for i in range(len(x_euler)-1):
u_euler[i+1] = u_euler[i] + dx*f(x_euler[i],u_euler[i])
plt.plot(x_ref,u_ref[:,0],'b',lw=3)
plt.plot(x_euler,u_euler[:,'-xr',lw=1)
处有 u_ref 的根，并且欧拉图距离不太远。