Use second-order and fourth-order Runge-Kutta methods

1 次查看（过去 30 天）
Cassidy Holene
Cassidy Holene 2021-3-6
回答： darova 2021-3-9
The one-dimensional linear convection equation for a scalar T is,
𝜕𝑇/𝜕𝑡 + 𝑎*𝜕𝑇/𝜕𝑥 = 0 0 ≤ 𝑥 ≤ 1
With the initial condition and boundary conditions of, (𝑥, 0) = 𝑒 ^(−200(𝑥−0.25)^ 2) , 𝑇(0,𝑡) = 0, 𝑇(1,𝑡): outflow condition Take 𝑎 = 1, ∆𝑥 = 0.01
Plot the results for t = 0.25, .5, 0.75 for both of them.
  2 个评论
darova
darova 2021-3-6
DO you have any attempts? Do you have difference scheme?
Cassidy Holene
Cassidy Holene 2021-3-6
a = 1; %initial condition
f = %given function
yprime = f(x,y) %direvative of f
xRange = [0,1]; %where the solution is sought on 0<=x<=1
yInitial = exp(-200*(xRange-0.25).^2); %column vector of initial values for y at x=0
numSteps = 100; %number of equally-sized steps to take from x1 to x2
t = 0.25; %first value of t changes to 0.5 and 0.75
x=zeros(numSteps+1,1); %initalzie x
x(1) = xRange(1); %set the first x value
h = ( xRange(2) - xRange(1) ) / numSteps; %time step
y(1,:) = (yInitial)';
%first try
for k = 1 : numSteps
x(k) = ??? ; %first x at k
y(k) = ??? ; %first y at k
x(k+1) = x(k) + h; %x at next k value
y(k+1,:) = y(k,:) + h * (feval( ??? ))'; %y at next k value
end
%second try
for i=1:(length(x)-1)
k1 = h*f(x(i),y(i)); %1st order
k2 = h*f(x(i)+0.5*k1,y(i)+0.5*h); %2nd order
k3 = h*f((x(i)+0.5*k2),(y(i)+0.5*h)); %3rd order
k4 = h*f((x(i)+k3),(y(i)+h)); %4th order
y(i+1) = y(i) + (1/6)*(k1+2*k2+2*k3+k4)*h; % main equation
end
figure, plot(x, y) % To see the solution results

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

darova
darova 2021-3-9
I'd try simple euler scheme first. The scheme i used
𝜕𝑇/𝜕𝑡 + 𝑎*𝜕𝑇/𝜕𝑥 = 0
clc,clear
a = 1;
x = 0:0.01:1;
t = 0:0.01:0.75;
U = zeros(length(t),length(x));
U(1,:) = exp(-200*(x-0.25).^2);
U(:,1) = 0;
dx = x(2)-x(1);
dt = t(2)-t(1);
for i = 1:size(U,1)-1
for j = 1:size(U,2)-1
U(i+1,j) = U(i,j) + a*dt/dx*(U(i,j+1)-U(i,j));
end
end
surf(U)

类别

在 Help Center 和 File Exchange 中查找有关 Vector Fields 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by