Clear Filters
Clear Filters

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

2 views (last 30 days)
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 Comments
darova
darova on 6 Mar 2021
DO you have any attempts? Do you have difference scheme?
Cassidy Holene
Cassidy Holene on 6 Mar 2021
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

Sign in to comment.

Answers (1)

darova
darova on 9 Mar 2021
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)

Categories

Find more on Vector Fields in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!