Wave Equation Using Matlab and Finite Differencing.

I'm trying to solve a wave equation, , over and ,with boundary/initial conditions:

I'm relatively new to numerically solving PDEs, but I'm not entirely sure what I'm doing wrong here. My code yields that the string vibrates to a height on the order of units, which is obviously unreasonable given the initial values. The spatial (n) / time (m) intervals are 25 and 700 respectively for the final equation, however they are shortened to 15 each in order to help me diagnose the code.

Here's what I have so far:

clear all;close all;clc;
% Initialize Variables
n = 15; % spacial steps
m = 15; % time steps
Li = 0; 
Lf = 2;
ti = 0;
tf = 50;
h = Lf/n; 
k = tf/m; 
u_left = 0;
u_right = 0;
u0 = @(x) -(x-1)^2 + 1;
ut0 = @(x) (x-1)^4 -1;
v = ut0;
c = sqrt(.25);
u = zeros(m+1,n+1);

% Boundary Conditions
u(:,1) = u_left;
u(:,n+1) = u_right;

% Initial Conditions
for i=2:n
    u(1,i) = u0(i*h); % Dirichlet Condition 
end
for i=2:n
    u(2,i) = u0(i*h) + ut0(i*h)*k; % Niemann Condition, used to fill timestep 2
end

% Solve using Finite Differencing
for i=3:m
    u(i,2:n) = c*k/h^2*u(i-1,1:n-1) + (1-2*c*k/h^2) * u(i-1,2:n) + c*k/h^2*u(i-1,3:n+1);
end
format long
u

u = 16×16

1.0e+31 * 0 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0 0 0 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 0 0 0 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 -0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 0 0 -0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 0 0 0.000000000000000 -0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0 0 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 0 0 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0 0 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 0 0 0.000000000000004 -0.000000000000005 0.000000000000004 -0.000000000000002 0.000000000000001 -0.000000000000000 0.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000001 0.000000000000002 -0.000000000000004 0.000000000000006 -0.000000000000004 0 0 -0.000000000001124 0.000000000001568 -0.000000000001274 0.000000000000701 -0.000000000000266 0.000000000000067 -0.000000000000011 0.000000000000016 -0.000000000000094 0.000000000000353 -0.000000000000895 0.000000000001580 -0.000000000001903 0.000000000001347 0

surf(u)

Answer 1

[nick]

Hi Tyler,

I found some errors in the equations for the initial condition and the finite difference method that you used to solve the wave equation above.

Instead of (ck/h^2), it should be ((ck/h)^2) in the finite difference method. Furthermore, it seems that there are certain terms missing in the finite difference equation. In the initial conditions, when using u0 and ut0, the input argument should correspond to (i-1) instead of (i). Kindly refer to the attached code below for your reference:

clear all;close all;clc;
%%
% Initialize Variables
n = 25; % spacial steps
Li = 0; %inital
Lf = 2;  %final
h = Lf/n; %total number steps x
m = 700; % time steps
ti = 0; %initial
tf = 50; %final
k = tf/m; %total number steps time 
u_left = 0;
u_right = 0;
u0 = @(x) -(x-1)^2 + 1; %u(x,0)
ut0 = @(x) (x-1)^4 -1; %d/dt(u(t)) @ (x,0)
v = ut0;
c = sqrt(.25);
u = zeros(m+1,n+1); %u(t,x)
%%
% Boundary Conditions
u(:,1) = u_left; %u(0,t)  x(0)==u(:,1)
u(:,n+1) = u_right; %u(2,t) 2==n+1
%%
% Initial Conditions
for i=2:n
    %i*h misses the first coordinate
    u(1,i) = u0((i-1)*h); % Dirichlet Condition 
end
for i=2:n
    %i*h misses the first coordinate
    u(2,i) = u0((i-1)*h) + ut0((i-1)*h)*k; % Niemann Condition, used to fill timestep 2
end
% Solve using Finite Differencing
r = (c*k/h)^2;
for i=3:m
    %u(i,2:n) = c*k/h^2*u(i-1,1:n-1) + (1-2*c*k/h^2) * u(i-1,2:n) + c*k/h^2*u(i-1,3:n+1);
    u(i,2:n) = 2*u(i-1,2:n) - u(i-2,2:n) + r*(u(i-1,3:n+1) - 2*u(i-1,2:n) + u(i-1,1:n-1));
end
format long
u

u = 701x26

0 0.153600000000000 0.294400000000000 0.422400000000000 0.537600000000000 0.640000000000000 0.729600000000000 0.806400000000000 0.870400000000000 0.921600000000000 0.960000000000000 0.985600000000000 0.998400000000000 0.998400000000000 0.985600000000000 0.960000000000000 0.921600000000000 0.870400000000000 0.806400000000000 0.729600000000000 0.640000000000000 0.537600000000000 0.422400000000000 0.294400000000000 0.153600000000000 0 0 0.133342354285714 0.258533668571429 0.374801554285714 0.481443840000000 0.577828571428571 0.663394011428571 0.737648640000000 0.800171154285714 0.850610468571429 0.888685714285714 0.914186240000000 0.926971611428571 0.926971611428571 0.914186240000000 0.888685714285714 0.850610468571429 0.800171154285714 0.737648640000000 0.663394011428571 0.577828571428571 0.481443840000000 0.374801554285714 0.258533668571428 0.133342354285714 0 0 0.111460218775510 0.220888911486881 0.325284741224490 0.423243365131195 0.513500874635569 0.594933795451895 0.666559087580175 0.727534145306123 0.777156797201166 0.814865306122449 0.840238369212828 0.852995117900875 0.852995117900875 0.840238369212828 0.814865306122449 0.777156797201166 0.727534145306122 0.666559087580175 0.594933795451895 0.513500874635568 0.423243365131195 0.325284741224490 0.220888911486880 0.111460218775510 0 0 0.089173203230202 0.182241112515024 0.274485002894032 0.363508069949426 0.447414450824062 0.524518934089368 0.593346959745344 0.652634619220563 0.701328655372166 0.738586462485869 0.763776086275957 0.776476223885286 0.776476223885286 0.763776086275956 0.738586462485869 0.701328655372166 0.652634619220563 0.593346959745344 0.524518934089368 0.447414450824061 0.363508069949426 0.274485002894032 0.182241112515024 0.089173203230202 0 0 0.067662396640294 0.143429087836491 0.223043359405010 0.302753027043494 0.379972419230083 0.452454587393348 0.518233451533290 0.575623800221336 0.623221290600345 0.659902448384602 0.684824667859820 0.697426211883144 0.697426211883144 0.684824667859820 0.659902448384602 0.623221290600345 0.575623800221336 0.518233451533290 0.452454587393348 0.379972419230083 0.302753027043494 0.223043359405010 0.143429087836491 0.067662396640294 0 0 0.047766763550846 0.105383880037009 0.171620728206717 0.241501676051645 0.311586266139109 0.379054282465639 0.441448125031237 0.496661279550187 0.542938317451062 0.578874895876718 0.603417757684298 0.615864731445231 0.615864731445231 0.603417757684298 0.578874895876718 0.542938317451062 0.496661279550187 0.441448125031237 0.379054282465639 0.311586266139109 0.241501676051645 0.171620728206717 0.105383880037009 0.047766763550846 0 0 0.029834290724335 0.069056571568591 0.120924360495992 0.180290910647037 0.242678633902541 0.304642702473969 0.363231698392337 0.415917804514786 0.460594506555602 0.495576593086214 0.519600155535194 0.531822588188255 0.531822588188255 0.519600155535194 0.495576593086214 0.460594506555602 0.415917804514785 0.363231698392337 0.304642702473969 0.242678633902541 0.180290910647037 0.120924360495992 0.069056571568591 0.029834290724335 0 0 0.013772829959352 0.035249493505774 0.071722484419444 0.119682260417902 0.173686567935887 0.229558475668472 0.283838834852145 0.333578067505206 0.376318573627655 0.410094273219492 0.433430606280716 0.445344532811329 0.445344532811329 0.433430606280717 0.410094273219492 0.376318573627655 0.333578067505206 0.283838834852145 0.229558475668472 0.173686567935887 0.119682260417902 0.071722484419444 0.035249493505773 0.013772829959352 0 0 -0.000753268563310 0.004431160533687 0.024809907028480 0.060278276068776 0.105066711833432 0.154157055673220 0.203540931745106 0.249843495010412 0.290255968479800 0.322531403581841 0.344984588887963 0.356492050112453 0.356492050112453 0.344984588887963 0.322531403581841 0.290255968479800 0.249843495010413 0.203540931745106 0.154157055673220 0.105066711833432 0.060278276068776 0.024809907028479 0.004431160533687 -0.000753268563310 0 0 -0.014095993030586 -0.023358968237627 -0.019095331685527 0.002731766752406 0.037304219425797 0.078814136202541 0.122628927712011 0.164935036634736 0.202571664039983 0.233010974593919 0.254357105439400 0.265346148004998 0.265346148004998 0.254357105439400 0.233010974593919 0.202571664039983 0.164935036634736 0.122628927712011 0.078814136202541 0.037304219425797 0.002731766752406 -0.019095331685527 -0.023358968237628 -0.014095993030586 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

surf(u)

I hope this helps.