Need help coding Finite Difference Scheme for a partial differential equation with both time and space
16 views (last 30 days)
Show older comments
Hi I need help coding a finite difference scheme for a partial differential equation. My math is correct and I have checked it with my professor but I cannot for the life of me get this code to work . I've included a picture of what the graph is supposed to look like. Also Included my code below.
Question:
Use a finite difference and Runge Kutta scheme to solve the following equation and then code the scheme and graph t=0 and t=2
Equation
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793349/image.png)
Derived Finite difference scheme (Checked with Professor)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793354/image.png)
Runge Kutta Scheme for Time
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793359/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793364/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793369/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793374/image.png)
Note: n is meant to represent the time space not an exponential
f represents the spatial discretization function
---
Solition Equation
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793379/image.png)
---
Initial Conditions
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793384/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793389/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793394/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793399/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793404/image.png)
Code:
clear; clc; close all;
global Nx Nt dx
%% Initial variables
tspan = [0 2];
dx = 0.1;
x = -10:dx:10;
Nx = length(x);
v = 20;
x0 = 0;
CFL = (4*sqrt(2)*sqrt(3))/9;
fprintf('%d\n',CFL)
dt = dx^3*CFL;
dt = 0.001;
t = 0:dt:2;
Nt = length(t);
u = zeros(length(x),length(t));
u1 = u;
%% solition solution
for i = 1:length(t)
for j = 1:length(x)
u1(j,i) = -v./(2*cosh(1/2.*sqrt(v).*(x(j)-v.*t(i)-x0))^2);
end
end
%% initial values case 1
u0 = u1(:,1); %u1(x,0);
u(:,1) = u0;
Loop for Runge Kutta Time space
for j = 1:Nt-1
uj1 = u(:,j);
a1 = dt*f5(uj1);
a2 = dt*f5(uj1+1/2*a1);
a3 = dt*f5(uj1+1/2*a2);
a4 = dt*f5(uj1+a3);
u(:,j+1) = uj1+1/6*(a1+2*(a2+a3)+a4);
%peridoic boundary
u(1,j+1) = u0(1);
u(end,j+1) = u0(end);
end
%% Graph
figure()
hold on
plot(x,u1(:,1))
plot(x,u(:,10),'r')
hold off
Function For Spacial Disc
I used an iterative function but I tried other ones too that I will attach below
function fufu = f5(u)
global Nx dx
u2 = zeros(Nx,1);
for i = 1:Nx
if (i == 1) %use forward diff
a1 = (u(i+1)-u(i))/dx;
b1 = 1/(2*dx^3)*(-5*u(i)+18*u(i+1)-24*u(i+2)+14*u(i+3)-3*u(i+4));
u2(i) = 6*u(i)*a1 - b1;
elseif i == 2
a1 = (u(i+1)-u(i-1))/(2*dx);
b1 = 1/(2*dx^3)*(-5*u(i)+18*u(i+1)-24*u(i+2)+14*u(i+3)-3*u(i+4));
u2(i) = 6*u(i)*a1 - b1;
elseif (i == Nx)
a1 = (-u(i-1)+u(i))/dx;
b1 = 1/(2*dx^3)*(5*u(i)-18*u(i-1)+24*u(i-2)-14*u(i-3)+3*u(i-4));
u2(i) = 6*u(i)*a1 - b1;
elseif i == Nx-1
a1 = (u(i+1)-u(i-1))/(2*dx);
b1 = 1/(2*dx^3)*(5*u(i)-18*u(i-1)+24*u(i-2)-14*u(i-3)+3*u(i-4));
u2(i) = 6*u(i)*a1 - b1;
else
a1 = (u(i+1)-u(i-1))/(2*dx);
b1 = (u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2))/(2*dx^3);
u2(i) = 6*u(i)*(a1) - b1;
end
end
Other Functions I've Tried:
I tried using matrix inversion, other forms of iterative and even a vectorized approach but nothing seems to work
function udot = f(u)
global Nx dx
u0 = u;
uj = zeros(size(u));
for i = 3:Nx-3
uj(i) = u(i)*6*((u(i+1)-u(i-1))/(2*dx)) - (u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2))/(2*dx^3);
end
uj(1:2) = u(1:2);
uj(Nx-1:Nx) = u(Nx-1:Nx);
uj(1) = uj(end);
udot = uj;
end
function uFun = f2(u)
global Nx dx
u0 = u;
a = 1/(2*dx);
dd(1:Nx) = 0;
du(1:Nx-1) = 1;
dl(1:Nx-1) = 1;
A = diag(dd)+diag(dl,-1)+diag(du,1);
A = sparse(A);
b = zeros(Nx,1);
b(1) = u0(end);
b(2:end) = u0(2:end)/(6*dx);
uFun = A\b;
end
%1st derivative
%6u*ux = whatever
function u2 = ufu(u)
global Nx dx
a = 1/(2*dx);
c = 1/dx^2;
dd(1:Nx) = 0;
d1(1:Nx-1) = dx^2*6*u(2:Nx) + 2;
d2(1:Nx-1) = -dx^2*6*u(2:Nx)- 2;
d3(1:Nx-2) = -1; %(u(i+1))
d4(1:Nx-2) = 1;
A = diag(dd)+diag(d1,1)+diag(d2,-1)+diag(d3,2)+diag(d4,-2);
A = sparse(A);
b = zeros(Nx,1);
b(1:Nx) = 2*dx^3*u(1:Nx);
u2 = A\b;
u2(1) = u2(end);
end
% Try just the simplified version
function ufun = f3(u)
global Nx dx
ut = diff(u); %du/dt
ut =[ut;0];
b = 1/(2*dx^3);
dd(1:Nx) = 0;
dL1(1:Nx-1) = -2*b;
dL2(1:Nx-2) = b;
dR1(1:Nx-1) = 2*b;
dR2(1:Nx-2) = -b;
A = diag(dd)+diag(dL1,-1)+diag(dL2,-2)+diag(dR1,1)+diag(dR2,2);
A = sparse(A);
A(:,1) = 0;
A(1,1) = 1;
A(:,end) = 0;
A(end,end) = 1;
b1 = zeros(size(u));
b1(1:Nx) = ut(1:Nx);
b1(1) = b(end);
ufun = A\b1;
end
% Trying iterative method again
function UT = f4(u)
global Nx dx
dudt = diff(u);
dudt = [dudt;0];
u2 = zeros(Nx,1);
for i = 1:length(u)
if i == 1
u2(i) = u(end);
% a1 = (u(i+1)-u(end-1))/(2*dx);
% u2(i) = 6*u(i)*a1 - (u(i+2)-2*u(i+1)+2*u(end-1)-u(end-2))/(2*dx^3);
elseif i == 2
a1 = (u(i+1)-u(i-1))/(2*dx);
u2(i) = 6*u(i)*a1 - (u(i+2)-2*u(i+1)+2*u(i-1))/(2*dx^3);
elseif i == Nx-1
a1 = (u(i+1)-u(i-1))/(2*dx);
u2(i) = 6*u(i)*a1 - (-2*u(i+1)+2*u(i-1)+u(i-2))/(2*dx^3);
elseif i == Nx
u2(i) = u(end);
% a1 = (u(2)-u(i-1))/(2*dx);
% u2(i) = 6*u(i)*a1 - (u(3)-2*u(2)+2*u(i-1)+u(i-2))/(2*dx^3);
elseif (i > 2) && (i < Nx-1)
a1 = (u(i+1)-u(i-1))/(2*dx);
b1 = (u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2))/(2*dx^3);
u2(i) = 6*u(i)*(a1) - b1;
end
end
% u2(1) = u2(end);
UT = u2;
end
What the correct code should look like
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793409/image.png)
0 Comments
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!