【Source Code】Scientific Computing MATLAB

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Scientific Computing MATLAB

Source code supporting the course "Introduction to Scientific and High Performance Computing with MATLAB" by Miguel O. Bernabeu

本帖隐藏的内容

• ScientificComputingMATLAB-master.zip (122.03 KB)
  2017-2-27 08:56:49 上传
• https://github.com/mobernabeu/ScientificComputingMATLAB
  

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏2 回帖

关键词：Scientific computing Comput MATLAB matla supporting course

1. function [ systemMatrix ] = AssemblyFiniteDifferencesMatrix( problemConfiguration )
   
2. %ASSEMBLYFINITEDIFFERENCESMATRIX Assembly a semi-implicit finite difference
   
3. %monodomain matrix (systemMatrix) according to the problem description in
   
4. %problemConfiguration.
   
5. %
   
6. %   systemMatrix is a tridiagonal matrix of the form
   
7. %
   
8. %     [a 2*b 0 ...       0]
   
9. %     [b  a  b  0 ...    0]
   
10. %     [0  b  a  b  0 ... 0]
    
11. %     [0 ... b  a  b ... 0]
    
12. %     [0 ... 0  b  a  b  0]
    
13. %     [0    ... 0  b  a  b]
    
14. %     [0       ... 0 2*b a]
    
15. %
    
16. %   with a = surfaceRatio*capacitance + 2*conductivity*timeStep/spaceStep^2,
    
17. %        b = -conductivity*timeStep/spaceStep^2
    
18. %
    
19. %        conductivity = { leftConductivity,  if columnIndex <= (numberCells+1)/2,
    
20. %                       { rightConductivity, otherwise.
    
21. 
    
22. % Introduction to Scientific and High Performance Computing with MATLAB
    
23. %
    
24. % Copyright (c) 2012-2014, Miguel O. Bernabeu, All rights reserved.
    
25. %
    
26. % This library is free software; you can redistribute it and/or
    
27. % modify it under the terms of the GNU Lesser General Public
    
28. % License as published by the Free Software Foundation; either
    
29. % version 3.0 of the License, or (at your option) any later version.
    
30. %
    
31. % This library is distributed in the hope that it will be useful,
    
32. % but WITHOUT ANY WARRANTY; without even the implied warranty of
    
33. % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    
34. % Lesser General Public License for more details.
    
35. %
    
36. % You should have received a copy of the GNU Lesser General Public
    
37. % License along with this library.
    
38. 
    
39. surfaceRatioTimesCapacitance = problemConfiguration.cellSurfaceAreaToVolumeRatio * problemConfiguration.cellMembraneCapacitance;
    
40. leftCondutivityTimesTimeStepOverSpaceStepSquared = problemConfiguration.leftConductivity * problemConfiguration.pdeTimeStep / problemConfiguration.spaceStep^2;
    
41. rightCondutivityTimesTimeStepOverSpaceStepSquared = problemConfiguration.rightConductivity * problemConfiguration.pdeTimeStep / problemConfiguration.spaceStep^2;
    
42. 
    
43. assert(problemConfiguration.numberCells > 1, 'Minimum problem size is currently 2');
    
44. 
    
45. systemMatrix = zeros(problemConfiguration.numberCells);
    
46. 
    
47. % First matrix row
    
48. systemMatrix(1,1) = surfaceRatioTimesCapacitance + 2*leftCondutivityTimesTimeStepOverSpaceStepSquared;
    
49. systemMatrix(1,2) = -2*leftCondutivityTimesTimeStepOverSpaceStepSquared;
    
50. 
    
51. % Rows 2 to numberCells-1
    
52. for rowIndex = 2:problemConfiguration.numberCells-1
    
53.     % Determine which conductivity to use (left half or right half)
    
54.     if rowIndex <= (problemConfiguration.numberCells+1)/2
    
55.         condutivityTimesTimeStepOverSpaceStepSquared = leftCondutivityTimesTimeStepOverSpaceStepSquared;
    
56.     else
    
57.         condutivityTimesTimeStepOverSpaceStepSquared = rightCondutivityTimesTimeStepOverSpaceStepSquared;
    
58.     end
    
59.    
    
60.     systemMatrix(rowIndex,rowIndex-1) = -condutivityTimesTimeStepOverSpaceStepSquared;
    
61.     systemMatrix(rowIndex,rowIndex) = surfaceRatioTimesCapacitance + 2*condutivityTimesTimeStepOverSpaceStepSquared;
    
62.     systemMatrix(rowIndex,rowIndex+1) = -condutivityTimesTimeStepOverSpaceStepSquared;
    
63. end
    
64. 
    
65. % Last matrix row
    
66. systemMatrix(problemConfiguration.numberCells,problemConfiguration.numberCells-1) = -2*rightCondutivityTimesTimeStepOverSpaceStepSquared;
    
67. systemMatrix(problemConfiguration.numberCells,problemConfiguration.numberCells) = surfaceRatioTimesCapacitance + 2*rightCondutivityTimesTimeStepOverSpaceStepSquared;
    
68. 
    
69. systemMatrix = sparse(systemMatrix);
    
70. 
    
71. end

复制代码

1. % Work out current time based on time step number
   
2. currentTime = problemConfiguration.pdeTimeStep*(timeStepNumber-1);
   
3. 
   
4. % Initialise stimuli vector to zero
   
5. stimuli = zeros(problemConfiguration.numberCells, 1);
   
6. 
   
7. % If its time to do it, stimulate two cells in the middle of the domain
   
8. if (currentTime>=problemConfiguration.stimulusStart && currentTime<problemConfiguration.stimulusStart+problemConfiguration.stimulusDuration)
   
9.     stimuli(floor(problemConfiguration.numberCells/2) : floor(problemConfiguration.numberCells/2) + 1) = problemConfiguration.stimulusMagnitude;
   
10. end
    
11. 
    
12. end

复制代码

1. function RegressionTest( problemConfiguration, conductionVelocities )
   
2. %REGRESSIONTEST This function checks the output of a simulation with 50
   
3. %cells and conductivities of 1.4 mS/cm against some known conduction
   
4. %velocity results
   
5. 
   
6. recordedConductionVelocity = 0.0521720;
   
7. 
   
8. if (problemConfiguration.numberCells==50 && problemConfiguration.leftConductivity==1.4 && problemConfiguration.rightConductivity==1.4)
   
9.    if any(conductionVelocities == Inf)
   
10.        warning('The simulation was not run for long enough to compute conduction velocities. Regression test could not to be performed.');
    
11.        return;
    
12.    end
    
13.    if any(abs(conductionVelocities - recordedConductionVelocity)>1e-4)
    
14.        error('Regression test failed. You may have introduced a bug while optimising the code! Review your last changes.');
    
15.    end
    
16. end
    
17. 
    
18. end

复制代码

kkkkkkkk

书么？

kankan

thanking you....