Date 14.5.99
Student 45834H, Mikko Hämäläinen

Problem Set 3

Problem 1

Problem 1 on Maple worksheet

Problem 2

Problem 2 on Maple worksheet

Problem 2.14

Files: main214.m, rhs214.m and exsol214.m



Here is the right hand side's function (rhs214.m):

function y=rhs214(x)
y=-exp(x-1);



Here is the function for exact solution (exsol214.m):

function y=exsol214(x)
y=exp(-1)*(exp(x)-1);



And here is the matlab code for the problem (main214.m):

clear;

n=5;
h=1/(n+1);

elem=((0:(n+1))')*h;

%% Scheme S1

maind=diag(ones(n+1,1));
upperd=diag(ones(n,1),-1);
lowerd=diag(ones(n,1),1);

A=(-lowerd+2*maind-upperd);
A(n+1,n+1)=1;

S1elem=((1:(n+1))')*h;
b=rhs214(S1elem)*h^2;
b(n+1)=h;

S1sol=A\b;
S1sol=[0; S1sol];

%% Scheme S2

mdiag=diag(ones(n+2,1));
updiag=diag(ones(n+1,1),1);
lwdiag=diag(ones(n+1,1),-1);

C=(-lwdiag+2*mdiag-updiag);
C(n+2,n+2)=1;
C(n+2,n+1)=0;
C(n+2,n)=-1;

S2elem=((1:(n+2))')*h;
d=rhs214(S2elem)*h^2;
d(n+2)=2*h;

S2sol=C\d;
S2sol=[0;S2sol(1:n+1)];

plot((0:n+1)'*h,S1sol,'g');
hold on;
plot((0:n+1)'*h,S2sol,'k');

N=15;
H=1/(N+1);
nodes=((0:(N+1))')*H;
solex=exsol214(nodes);
plot((0:N+1)'*H,solex,'ro');
hold off;

Problem 2.17

Files: ex217.m, f217.m and sol217.m



Here is the solution function of the problem (sol217.m):

function y=sol217(x)
y=1-(1-exp(-10))*x-exp(-10*x);



Here is the function of the problem (f217.m):

function y=f217(x)
y=100*exp(-10*x);


 
And here is the Matlab session of the problem (ex217.m):

clear;

for i=1:5

n=10*(2^(i-1))-1;
h=1/(n+1);

nodes=((0:(n+1)).')*h;

exactsol=sol217(nodes);

maind=diag(ones(n,1));
lowd=diag(ones(n-1,1),-1);
uppd=diag(ones(n-1,1),1);

A=(-lowd+2*maind-uppd);

discret=((1:(n)).')*h;
b=h^2*f217(discret);

solvect=A\b;
solution=[0;solvect;0];

H(i)=h;
error(i)=max(abs(solution-exactsol));
estimate(i)=10^4*h^2/96;

hold on;
if i==1
plot((0:n+1).'*h,exactsol,'r');
plot((0:n+1).'*h,solution,'g*');
end
if i==2
plot((0:n+1).'*h,solution,'bx');
end
if i==3
plot((0:n+1).'*h,solution,'y+');
end
if i==4
plot((0:n+1).'*h,solution,'c.');
end
if i==5
plot((0:n+1).'*h,solution,'k-.');
end
hold off;

end

H
error
estimate